17 List of examples

List of examples shipped in the distribution of Optimization Toolbox for MATLAB:

Table 17.1 List of distributed examples

File

Description

advs.m

Advanced simplex hot-start examples

affco1.m

A simple problem using affine conic constraints

affco2.m

A simple problem using affine conic constraints

callback.m

An example of data/progress callback

callback_handler.m

Log handler definition for callback.m

ceo1.m

A simple conic exponential problem

cqo1.m

A simple conic quadratic problem

feasrepairex1.m

A simple example of how to repair an infeasible problem

gp1.m

A simple geometric program (GP) in conic form

lo1.m

A simple linear problem using msklpopt

lo2.m

A simple linear problem using mosekopt

lo3.m

A simple linear problem using linprog

mico1.m

A simple mixed-integer conic problem

milo1.m

A simple mixed-integer linear problem

mioinitsol.m

A simple mixed-integer linear problem with an initial guess

normex.m

Demonstrates least squares and other norm minimization problems

opt_server_sync.m

Uses MOSEK OptServer to solve an optimization problem synchronously

parameters.m

Shows how to set optimizer parameters and read information items

portfolio_1_basic.m

Portfolio optimization - basic Markowitz model

portfolio_2_frontier.m

Portfolio optimization - efficient frontier

portfolio_3_impact.m

Portfolio optimization - market impact costs

portfolio_4_transcost.m

Portfolio optimization - transaction costs

portfolio_5_card.m

Portfolio optimization - cardinality constraints

pow1.m

A simple power cone problem

qcqo1.m

A simple quadratically constrained quadratic problem

qo1.m

A simple quadratic problem

qo2.m

A simple quadratic problem

reoptimization.m

Demonstrate how to modify and re-optimize a linear problem

response.m

Demonstrates proper response handling

rlo1.m

Robust linear optimization example, part 1

rlo2.m

Robust linear optimization example, part 2

sdo1.m

A simple semidefinite problem with one matrix variable and a quadratic cone

sdo2.m

A simple semidefinite problem with two matrix variables

sensitivity.m

Sensitivity analysis performed on a small linear problem

sensitivity2.m

Sensitivity analysis performed on a small linear problem

simple.m

A simple I/O example: read problem from a file, solve and write solutions

solutionquality.m

Demonstrates how to examine the quality of a solution

Additional examples can be found on the MOSEK website and in other MOSEK publications.

advs.m

Listing 17.1 advs.m Click here to download.
%%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      advs.m
%
%  Demonstrates hot-start capabilities for the simplex optimizers.
%%

clear prob param bas

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Initial problem (with hotstart) %%%%%%%
% 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)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Add a variable %%%%%%%%%%%%%%%%%
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'


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Fixing a variable %%%%%%%%%%%%%%%%%%
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'


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Add a constraint %%%%%%%%%%%%%%%%%%
% 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'


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% Remove a constraint %%%%%%%%%%%%%%%%%
% Modify the problem.
prob.a       = prob.a(1:end-1,:);
prob.blc     = prob.blc(1:end-1);
prob.buc     = prob.buc(1:end-1);

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

% Set the solution to the modified problem.
bas.skc      = bas.skc(1:end-1,:);
bas.xc       = bas.xc(1:end-1);
bas.y        = bas.y(1:end-1);
bas.slc      = bas.slc(1:end-1);
bas.suc      = bas.suc(1:end-1);

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

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

res.sol.bas.xx'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Remove a variable %%%%%%%%%%%%%%%%
% Modify the problem.
prob.c       = prob.c(1:end-1);
prob.a       = prob.a(:,1:end-1);
prob.blx     = prob.blx(1:end-1);
prob.bux     = prob.bux(1:end-1);

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

% Set the solution to the modified problem.
bas.xx       = bas.xx(1:end-1);
bas.skx      = bas.skx(1:end-1,:);
bas.slx      = bas.slx(1:end-1);
bas.sux      = bas.sux(1:end-1);

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

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

res.sol.bas.xx'

affco1.m

Listing 17.2 affco1.m Click here to download.
%%
%  File : affco1.m
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  Description :
%    Implements a basic tutorial example with affine conic constraints:
%
%    maximize x_1^(1/3) + (x_1+x_2+0.1)^(1/4)
%    st.      (x_1-0.5)^2 + (x_2-0.6)^2 <= 1
%             0 <= x_1 <= x_2 + 1
% 
%%
function affco1()

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

% Variables [x1; x2; t1; t2]
prob.c = [0, 0, 1, 1];

% Linear inequality x_1 - x_2 <= 1
prob.a = sparse([1, -1, 0, 0]);
prob.buc = 1;
prob.blc = [];

% The quadratic cone
FQ = sparse([zeros(1,4); speye(2) zeros(2,2)]);
gQ = [1 -0.5 -0.6]';
cQ = [res.symbcon.MSK_CT_QUAD 3];

% The power cone for (x_1, 1, t_1) \in POW3^(1/3,2/3)
FP1 = sparse([1 0 zeros(1,2); zeros(1,4); zeros(1,2) 1 0]);
gP1 = [0 1 0]';
cP1 = [res.symbcon.MSK_CT_PPOW 3 2 1/3 2/3];

% The power cone for (x_1+x_2+0.1, 1, t_2) \in POW3^(1/4,3/4)
FP2 = sparse([1 1 zeros(1,2); zeros(1,4); zeros(1,2) 0 1]);
gP2 = [0.1 1 0]';
cP2 = [res.symbcon.MSK_CT_PPOW 3 2 1.0 3.0];

% All cones
prob.f = [FQ; FP1; FP2];
prob.g = [gQ; gP1; gP2];
prob.cones = [cQ cP1 cP2];

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

res.sol.itr.pobjval
res.sol.itr.xx(1:2)

affco2.m

Listing 17.3 affco2.m Click here to download.
%%
%  File : affco2.m
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  Description :
%    Implements a basic tutorial example with affine conic constraints:
%
%    minimize t
%    st.      (d, z1*y1,... zn*yn) \in Q^{n+1}
%             (yi, 1, ai*t)        \in EXP, i=1,\ldots,n
%
%    with input ai<0, zi, d.
%
%    See also https://docs.mosek.com/modeling-cookbook/expo.html#hitting-time-of-a-linear-system
% 
%%
function affco2()

n = 2;
z = [2.2 1.3]';
a = [-0.3 -0.06]';
d = 0.5;

t = firstHittingTime(n, z, a, d);
disp(sprintf('\nt = %.4e', t))


function t = firstHittingTime(n, z, a, d)

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

% Variables [t, y1, ..., yn]
prob.a = sparse(0, n+1);
prob.c = [1 zeros(1,n)];

% Quadratic cone
FQ = diag([0; z]);
gQ = [d; zeros(n,1)];

% All exponential cones
FE = sparse([1:3:3*n    3:3:3*n], ...
            [2:n+1      ones(1,n)], ...
            [ones(1,n)  a']);
gE = repmat([0; 1; 0], n, 1);

% Assemble input data
prob.f = [FQ; FE];
prob.g = [gQ; gE];
prob.cones = [res.symbcon.MSK_CT_QUAD n+1 repmat([res.symbcon.MSK_CT_PEXP 3], 1, n)];

% Solve
[r, res] = mosekopt('minimize', prob);
t = res.sol.itr.xx(1)

callback.m

Listing 17.4 callback.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      callback.m
%
%  Purpose :   To demonstrate how to use the progress
%              callback.

% Example of how to attach the callback function
function callback(arg1, arg2)

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 = []; 

% Define user defined handle.
[r,res]             = mosekopt('echo(0) symbcon');                 
data.maxtime        = 100.0;
data.symbcon        = res.symbcon;

callback.iter       = 'callback_handler';	% Defined in callback_handler.m
callback.iterhandle = data;

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

callback_handler.m

Listing 17.5 callback_handler.m Click here to download.
% The callback function. It will display some information
% during optimization.
% 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.
function [r] = callback_handler(handle,where,info)

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

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

if handle.symbcon.MSK_CALLBACK_END_INTPNT==where
  disp(sprintf('Interior-point optimizer terminated\n'));
  disp(sprintf('Interior-point primal obj.: %e\n', info.MSK_DINF_INTPNT_PRIMAL_OBJ));
  disp(sprintf('Iterations: %d\n', info.MSK_IINF_INTPNT_ITER));
end

if handle.symbcon.MSK_CALLBACK_NEW_INT_MIO==where
  disp(sprintf('New mixed-integer solution found\n'));
  disp(sprintf('Best objective.: %e\n', info.MSK_DINF_MIO_OBJ_BOUND));
end

% Decide if to terminate the optimization
% Terminate when cputime > handle.maxtime
if info.MSK_DINF_INTPNT_TIME > handle.maxtime
  r = 1;
else
  r = 0;
end

ceo1.m

Listing 17.6 ceo1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      ceo1.m
%
%  Purpose:   To demonstrate how to solve a small conic exponential
%              optimization problem using the MOSEK Matlab Toolbox.
%

function ceo1()

clear prob;

[r, res] = mosekopt('symbcon');
% Specify the non-conic part of the problem.

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

% Specify the cones.

prob.cones.type   = [res.symbcon.MSK_CT_PEXP];
prob.cones.sub    = [1, 2, 3];
prob.cones.subptr = [1];
% The field 'type' specifies the cone types, in this case an exponential
% cone with key MSK_CT_PEXP.
%
% The fields 'sub' and 'subptr' specify the members of the cones,
% i.e., the above definitions imply that 
%   x(1) >= x(2)*exp(x(3)/x(2))

% Optimize the problem. 

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

% Display the primal solution.

res.sol.itr.xx'

cqo1.m

Listing 17.7 cqo1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      cqo1.m
%
%  Purpose:   To demonstrate how to solve a small conic quadratic
%              optimization problem using the MOSEK Matlab Toolbox.
%

function cqo1()

clear prob;

[r, res] = mosekopt('symbcon');
% Specify the non-conic part of the problem.

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

% Specify the cones.

prob.cones.type   = [res.symbcon.MSK_CT_QUAD, res.symbcon.MSK_CT_RQUAD];
prob.cones.sub    = [4, 1, 2, 5, 6, 3];
prob.cones.subptr = [1, 4];
% The field 'type' specifies the cone types, i.e., quadratic cone
% or rotated quadratic cone. The keys for the two cone types are MSK_CT_QUAD
% and MSK_CT_RQUAD, respectively.
%
% The fields 'sub' and 'subptr' specify the members of the cones,
% i.e., the above definitions imply that 
%   x(4) >= sqrt(x(1)^2+x(2)^2) and 2 * x(5) * x(6) >= x(3)^2.

% Optimize the problem. 

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

% Display the primal solution.

res.sol.itr.xx'

feasrepairex1.m

Listing 17.8 feasrepairex1.m Click here to download.
%%
%  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File:      feasrepairex1.m
%
%  Purpose:   To demonstrate how to use the MSK_primalrepair function to
%             repair an infeasible problem.
%
%%

function feasrepairex1(inputfile)

cmd = sprintf('read(%s)', inputfile);
[r,res]=mosekopt(cmd);

res.prob.primalrepair = [];
res.prob.primalrepair.wux = [1,1];
res.prob.primalrepair.wlx = [1,1];
res.prob.primalrepair.wuc = [1,1,1,1];
res.prob.primalrepair.wlc = [1,1,1,1];

param.MSK_IPAR_LOG_FEAS_REPAIR = 3;
[r,res]=mosekopt('minimize primalrepair',res.prob,param);
fprintf('Return code: %d\n',r);

end

gp1.m

Listing 17.9 gp1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      gp1.m
%
%  Purpose:    Demonstrates how to solve a simple Geometric Program (GP)
%              cast into conic form with exponential cones and log-sum-exp.
%
%              Example from
%                https://gpkit.readthedocs.io/en/latest/examples.html//maximizing-the-volume-of-a-box
%

%
% maximize     h*w*d
% subjecto to  2*(h*w + h*d) <= Awall
%                       w*d <= Afloor
%              alpha <= h/w <= beta
%              gamma <= d/w <= delta
%
% Variable substitutions:  h = exp(x), w = exp(y), d = exp(z).
%
% maximize     x+y+z
% subject      log( exp(x+y+log(2/Awall)) + exp(x+z+log(2/Awall)) ) <= 0
%                              y+z <= log(Afloor)  
%              log( alpha ) <= x-y <= log( beta )
%              log( gamma ) <= z-y <= log( delta )
%
%
% log(sum(exp(A*x+b))) <= 0 is equivalent to sum(u) == 1, (ui,1,ai*x+bi) in Kexp,
% so we have two exp-cones and two auxilliary variables u1,u2.
%
% We order variables as (x,y,z,u1,u2),

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

% Input data
Awall = 200;
Afloor = 50;
alpha = 2;
beta = 10;
gamma = 2;
delta = 10;

% Objective
prob = [];
prob.c = [1, 1, 1, 0, 0]';

% Linear constraints:
% [ 0  0  0  1  1 ]                    == 1
% [ 0  1  1  0  0 ]                    <= log(Afloor)     
% [ 1 -1  0  0  0 ]                    in [log(alpha), log(beta)]
% [ 0 -1  1  0  0 ]                    in [log(gamma), log(delta)]
%
prob.a = [ 0  0  0  1  1;
           0  1  1  0  0;
           1 -1  0  0  0;
           0 -1  1  0  0 ];
           
prob.blc = [ 1; -inf;        log(alpha); log(gamma) ];
prob.buc = [ 1; log(Afloor); log(beta);  log(delta) ];

prob.blx = [ -inf; -inf; -inf; -inf; -inf];
prob.bux = [ inf; inf; inf; inf; inf];

% The conic part FX+g \in Kexp x Kexp
%   x  y  z  u  v
% [ 0  0  0  1  0 ]    0 
% [ 0  0  0  0  0 ]    1               in Kexp
% [ 1  1  0  0  0 ]    log(2/Awall)
%
% [ 0  0  0  0  1 ]    0
% [ 0  0  0  0  0 ]    1               in Kexp
% [ 1  0  1  0  0 ] +  log(2/Awall)
%
%         
prob.f = sparse([0 0 0 1 0; 
                 0 0 0 0 0;
                 1 1 0 0 0;
                 0 0 0 0 1;
                 0 0 0 0 0;
                 1 0 1 0 0]);
            
prob.g = [ 0; 1; log(2/Awall); 0; 1; log(2/Awall)];

prob.cones = [ res.symbcon.MSK_CT_PEXP, 3, res.symbcon.MSK_CT_PEXP, 3 ];

% Optimize and print results
[r,res]=mosekopt('maximize',prob);
exp(res.sol.itr.xx(1:3))

lo1.m

Listing 17.10 lo1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      lo1.m
%
%  Purpose : Demonstrates a simple linear problem.
%
function lo1()

c     = [3 1 5 1]';
a     = [[3 1 2 0];[2 1 3 1];[0 2 0 3]];
blc   = [30 15  -inf]';
buc   = [30 inf 25 ]';
blx   = zeros(4,1);
bux   = [inf 10 inf inf]';

[res] = msklpopt(c,a,blc,buc,blx,bux,[],'maximize');
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.

lo2.m

Listing 17.11 lo2.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      lo2.m
%
%  Purpose : Demonstrates a simple linear problem.
%

function lo2()
clear prob;

% Specify the c vector.
prob.c  = [3 1 5 1]';

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

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

% Specify lower bounds of the constraints.
prob.blc = [30 15  -inf]';

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

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

% Specify upper bounds of the variables.
prob.bux = [inf 10 inf inf]';

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

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

lo3.m

Listing 17.12 lo3.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      lo3.m
%
%  Purpose : Demonstrates a simple linear problem solved with MOSEK's linprog.
%
f     = - [3 1 5 1]';				% minus because we maximize
A     = [[-2 -1 -3 -1]; [0 2 0 3]];
b     = [-15  25]';
Aeq   = [3 1 2 0];
beq   = 30;
l     = zeros(4,1);
u     = [inf 10 inf inf]';

% Example of setting options for linprog
% Get default options
opt = mskoptimset('');
% Turn on diagnostic output
opt = mskoptimset(opt,'Diagnostics','on');
% Set a MOSEK option, in this case turn basic identification off.
opt = mskoptimset(opt,'MSK_IPAR_INTPNT_BASIS','MSK_OFF');
% Modify a MOSEK parameter with double value
opt = mskoptimset(opt,'MSK_DPAR_INTPNT_TOL_INFEAS',1e-12);

[x,fval,exitflag,output,lambda] = linprog(f,A,b,Aeq,beq,l,u,opt);

x
fval
exitflag
output
lambda

mico1.m

Listing 17.13 mico1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File      : mico1.m
%
%  Purpose   : Demonstrates how to solve a small mixed
%              integer conic optimization problem.
%
%              minimize    x^2 + y^2
%              subject to  x >= e^y + 3.8
%                          x, y - integer
%

function mico1()

[rcode, res] = mosekopt('symbcon echo(0)');
symbcon = res.symbcon;
clear prob       

% The full variable is [t; x; y]
prob.c        = [1 0 0];
prob.a        = sparse(0,3);   % No constraints

% Conic part of the problem
prob.f = sparse([ eye(3); 
                  0 1 0; 
                  0 0 0; 
                  0 0 1 ]);
prob.g = [0 0 0 -3.8 1 0]';
prob.cones = [symbcon.MSK_CT_QUAD 3 symbcon.MSK_CT_PEXP 3];

% Specify indexes of variables that are integers
prob.ints.sub = [2 3];

% It is as always possible (but not required) to input an initial solution
% to start the mixed-integer solver. 
prob.sol.int.xx = [0, 9, -1];

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

% The integer solution (x,y) 
res.sol.int.xx(2:3)

milo1.m

Listing 17.14 milo1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File      : milo1.m
%
%  Purpose   : Demonstrates how to solve a small mixed
%              integer linear optimization problem.
%

function milo1()
clear prob       
prob.c        = [1 0.64];
prob.a        = [[50 31];[3 -2]];
prob.blc      = [-inf -4];
prob.buc      = [250 inf];
prob.blx      = [0 0];
prob.bux      = [inf inf];

% Specify indexes of variables that are integer
% constrained.

prob.ints.sub = [1 2];

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

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

mioinitsol.m

Listing 17.15 mioinitsol.m Click here to download.
%%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      mioinitsol.m
%
%   Purpose:   To demonstrate how to solve a MIP with a start
%   guess.
%%

function mioinitsol()
[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];

% Specify start guess for the integer variables.
prob.sol.int.xx  = [1 1 0 nan]';

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

try
  % Display the optimal solution.
  res.sol.int.xx'

  % Was the initial solution used ?
  res.info.MSK_IINF_MIO_CONSTRUCT_SOLUTION
  res.info.MSK_DINF_MIO_CONSTRUCT_SOLUTION_OBJ

catch
  fprintf('MSKERROR: Could not get solution')
end

normex.m

Listing 17.16 normex.m Click here to download.
%%
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      normex.m
%
%  Purpose:   Demonstrates various norm minimization problems
%
%             * least squares regression
%             * ridge regularization
%             * lasso regularization
%             * p-norm minimization
%
%%

% This function runs a few examples
function normex()
% Create some random data
N  = 100;
x1 = sort(rand(N,1));
p  = @(t) 1.5*t.^3-2*t.^2+0.5*t.^1-1;
y1 = p(x1) + 0.01*randn(N,1);
xp = smallVander(x1, 5);
figure(1);
scatter(x1, y1);
hold on

% Least squares regression
x = norm_lse(xp, y1, zeros(0,6), []);
plot(0:0.01:1, smallVander((0:0.01:1)', 5)*x);

% With ridge regularization
x = norm_lse_reg(xp, y1, zeros(0,6), [], 0.1);
plot(0:0.01:1, smallVander((0:0.01:1)', 5)*x);

% Quadratic ridge
x = norm_lse_reg_quad(xp, y1, zeros(0,6), [], 0.2);
plot(0:0.01:1, smallVander((0:0.01:1)', 5)*x);

% Completely random large data for lasso example
N = 100;
K = 3000;
F = sprandn(K, N, 0.5);
g = randn(K, 1);
disp(sprintf('Lasso regularization'));
for gamma=[0.01, 0.1, 0.3, 0.6, 0.9, 1.3]
   x = norm_lse_lasso(F, g, zeros(0,N), [], gamma);
   disp(sprintf('Gamma %.4f  density %.0f%%   |Fx-g|_2: %.4f', gamma, sum(abs(x)>1e-6)/N*100, norm(F*x-g)));
end

% Example with the p-norm cone for various p
% We add a far outlier to the first example
x12 = [x1; 0.73];
y12 = [y1; -0.99];
xp2 = [x12.^5, x12.^4, x12.^3, x12.^2, x12.^1, x12.^0];
figure(2);
parts = [];
labels = [];
scatter(x12, y12);
hold on
for p=[1.1, 2.0, 3.0, 6.0]
   x = norm_p_norm(xp2, y12, zeros(0,6), [], p);
   parts = [parts; plot(0:0.01:1, smallVander((0:0.01:1)', 5)*x);];
   labels = [labels; sprintf('%.1f', p)];
end
legend(parts, labels);
end

% Least squares regression
% minimize \|Fx-g\|_2
function x = norm_lse(F,g,A,b)
clear prob;
[r, res] = mosekopt('symbcon');
n = size(F,2);
k = size(g,1);
m = size(A,1);

% Linear constraints in [x; t]
prob.a   = [A, zeros(m,1)];
prob.buc = b;
prob.blc = b;
prob.blx = -inf*ones(n+1,1);
prob.bux = inf*ones(n+1,1);
prob.c   = [zeros(n,1); 1];

% Affine conic constraint
prob.f = sparse([zeros(1,n), 1; F, zeros(k,1)]);
prob.g = [0; -g];
prob.cones = [ res.symbcon.MSK_CT_QUAD k+1 ];

% Solve
[r, res] = mosekopt('minimize echo(0)', prob);
x = res.sol.itr.xx(1:n);
end

% Least squares regression with regularization
% minimize \|Fx-g\|_2 + gamma*\|x\|_2
function x = norm_lse_reg(F,g,A,b,gamma)
clear prob;
[r, res] = mosekopt('symbcon');
n = size(F,2);
k = size(g,1);
m = size(A,1);

% Linear constraints in [x; t1; t2]
prob.a   = [A, zeros(m,2)];
prob.buc = b;
prob.blc = b;
prob.blx = -inf*ones(n+2,1);
prob.bux = inf*ones(n+2,1);
prob.c   = [zeros(n,1); 1; gamma];

% Affine conic constraint
prob.f = sparse([zeros(1,n),        1, 0; ...
                 F,           zeros(k,2); ...
                 zeros(1,n),        0, 1; ...
                 eye(n),      zeros(n,2) ]);
prob.g = [0; -g; zeros(n+1,1)];
prob.cones = [ res.symbcon.MSK_CT_QUAD k+1 res.symbcon.MSK_CT_QUAD n+1 ];

% Solve
[r, res] = mosekopt('minimize echo(0)', prob);
x = res.sol.itr.xx(1:n);
end

% Least squares regression with regularization
% The "classical" quadratic version
% minimize \|Fx-g\|_2^2 + gamma*\|x\|_2^2
function x = norm_lse_reg_quad(F,g,A,b,gamma)
clear prob;
[r, res] = mosekopt('symbcon');
n = size(F,2);
k = size(g,1);
m = size(A,1);

% Linear constraints in [x; t1; t2]
prob.a   = [A, zeros(m,2)];
prob.buc = b;
prob.blc = b;
prob.blx = -inf*ones(n+2,1);
prob.bux = inf*ones(n+2,1);
prob.c   = [zeros(n,1); 1; gamma];

% Affine conic constraint
prob.f = sparse([zeros(1,n),        1, 0; ...
                 zeros(1,n+2)           ; ...
                 F,           zeros(k,2); ...
                 zeros(1,n),        0, 1; ...
                 zeros(1,n+2)           ; ...                
                 eye(n),      zeros(n,2) ]);
prob.g = [0; 0.5; -g; 0; 0.5; zeros(n,1)];
prob.cones = [ res.symbcon.MSK_CT_RQUAD k+2 res.symbcon.MSK_CT_RQUAD n+2 ];

% Solve
[r, res] = mosekopt('minimize echo(0)', prob);
x = res.sol.itr.xx(1:n);
end

% Least squares regression with lasso regularization
% minimize \|Fx-g\|_2 + gamma*\|x\|_1
function x = norm_lse_lasso(F,g,A,b,gamma)
clear prob;
[r, res] = mosekopt('symbcon');
n = size(F,2);
k = size(g,1);
m = size(A,1);

% Linear constraints in [x; u; t1; t2]
prob.a   = [A,         zeros(m,n+2)      ; ...
            eye(n),    eye(n), zeros(n,2); ...
            -eye(n),   eye(n), zeros(n,2); ...
            zeros(1,n) -ones(1,n), 0, 1 ];
prob.buc = [b; inf*ones(2*n+1,1)];
prob.blc = [b; zeros(2*n+1,1)];
prob.blx = -inf*ones(2*n+2,1);
prob.bux = inf*ones(2*n+2,1);
prob.c   = [zeros(2*n,1); 1; gamma];

% Affine conic constraint
prob.f = sparse([zeros(1,2*n), 1, 0; F, zeros(k,n+2)]);
prob.g = [0; -g];
prob.cones = [ res.symbcon.MSK_CT_QUAD k+1 ];

% Solve
[r, res] = mosekopt('minimize echo(0)', prob);
x = res.sol.itr.xx(1:n);
end

% P-norm minimization
% minimize \|Fx-g\|_p
function x = norm_p_norm(F,g,A,b,p)
clear prob;
[r, res] = mosekopt('symbcon');
n = size(F,2);
k = size(g,1);
m = size(A,1);

% Linear constraints in [x; r; t]
prob.a   = [A, zeros(m,k+1); zeros(1,n), ones(1,k), -1];
prob.buc = [b; 0];
prob.blc = [b; 0];
prob.blx = -inf*ones(n+k+1,1);
prob.bux = inf*ones(n+k+1,1);
prob.c   = [zeros(n+k,1); 1];

% Permutation matrix which picks triples (r_i, t, F_ix-g_i)
M = [];
for i=1:3
  M = [M, sparse(i:3:3*k, 1:k, ones(k,1), 3*k, k)];
end

% Affine conic constraint
prob.f = M * sparse([zeros(k,n), eye(k), zeros(k,1); zeros(k,n+k), ones(k,1); F, zeros(k,k+1)]);
prob.g = M * [zeros(2*k,1); -g];
prob.cones = [ repmat([res.symbcon.MSK_CT_PPOW, 3, 2, 1.0, p-1], 1, k) ];

% Solve
[r, res] = mosekopt('minimize echo(0)', prob);
x = res.sol.itr.xx(1:n);
end

% Just a few columns of the Vandermode matrix
function M = smallVander(v, k)
M = [];
for i=0:k 
    M = [v.^i, M];
end
end

opt_server_sync.m

Listing 17.17 opt_server_sync.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      opt_server_sync.m
%
%  Purpose : Demonstrates solving a problem remotely using the OptServer.
%
%

function opt_server_sync(inputfile, addr)
clear prob;
clear param;
clear optserver;

% We read some problem from a file or set it up here
cmd = sprintf('read(%s)', inputfile);
[r,res] = mosekopt(cmd);
prob = res.prob;

% OptServer data
% "host" should be 'http://server:port`
optserver.host = addr;

% Perform the optimization with full log output.
[r,res] = mosekopt(sprintf('%s echo(10)', prob.objsense), prob, [], [], optserver); 

% Use the optimal x solution.
xx = res.sol.bas.xx;

parameters.m

Listing 17.18 parameters.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      parameters.m
%
%  Purpose:    Demonstrates a very simple example about how to get/set
%              parameters with MOSEK Matlab Toolbox
%

function r = parameters()

fprintf('Test MOSEK parameter get/set functions\n\n');

% Set up a small problem
prob.c  = [3 1 5 1]';
subi   = [1 1 1 2 2 2 2 3 3];
subj   = [1 2 3 1 2 3 4 2 4];
valij  = [3 1 2 2 1 3 1 2 3];
prob.a = sparse(subi,subj,valij);
prob.blc = [30 15  -inf]';
prob.buc = [30 inf 25 ]';
prob.blx = zeros(4,1);
prob.bux = [inf 10 inf inf]';

% Get default parameter values
[r,resp]=mosekopt('param'); 
fprintf('Default value for parameter MSK_IPAR_LOG = %d\n', resp.param.MSK_IPAR_LOG);
    
% Set log level (integer parameter)
param.MSK_IPAR_LOG = 1;
% Select interior-point optimizer... (integer parameter)
param.MSK_IPAR_OPTIMIZER = 'MSK_OPTIMIZER_INTPNT';
% ... without basis identification (integer parameter)
param.MSK_IPAR_INTPNT_BASIS = 'MSK_BI_NEVER';
% Set relative gap tolerance (double parameter)
param.MSK_DPAR_INTPNT_CO_TOL_REL_GAP = 1.0e-7;

% Use in mosekopt
[r,resp] = mosekopt('minimize', prob, param);

% Demonstrate information items after optimization
[r,res] = mosekopt('minimize info', prob);

res.info.MSK_DINF_OPTIMIZER_TIME
res.info.MSK_IINF_INTPNT_ITER

end

portfolio_1_basic.m

Listing 17.19 portfolio_1_basic.m Click here to download.
%%
%  File : portfolio_1_basic.m
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  Description :
%    Implements a basic portfolio optimization model.
% 
%%
function portfolio_1_basic()

n      = 3;
w      = 1.0;   
mu     = [0.1073,0.0737,0.0627]';
x0     = [0.0,0.0,0.0]';
gammas = [0.035,0.040,0.050,0.060,0.070,0.080 ,0.090]';
GT     = [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238  ; ...
           0.0              , 0.102863378954911  , -0.00222873156550421 ; ...
           0.0              , 0.0                ,  0.0338148677744977 ];
         
disp('Basic Markowitz portfolio optimization')
for gamma = gammas'
    er = BasicMarkowitz(n,mu,GT,x0,w,gamma);
    disp(sprintf('Expected return: %.4e Std. deviation: %.4e', er, gamma));
end


%
%    Purpose:
%        Computes the optimal portfolio for a given risk 
%     
%    Input:
%        n: Number of assets
%        mu: An n dimmensional vector of expected returns
%        GT: A matrix with n columns so (GT')*GT  = covariance matrix
%        x0: Initial holdings 
%        w: Initial cash holding
%        gamma: Maximum risk (=std. dev) accepted
%     
%    Output:
%        Optimal expected return and the optimal portfolio     
%
function er = BasicMarkowitz(n,mu,GT,x0,w,gamma)

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

% Objective vector - expected return
prob.c = mu;

% The budget constraint  - e'x = w + sum(x0)
prob.a = ones(1,n);
prob.blc = w + sum(x0);
prob.buc = w + sum(x0);

% Bounds exclude shortselling
prob.blx = zeros(n,1);
prob.bux = inf*ones(n,1); 

% An affine conic constraint: [gamma, GT*x] in quadratic cone
prob.f = sparse([ zeros(1,n); GT ]);
prob.g = [gamma; zeros(n,1)];
prob.cones = [ res.symbcon.MSK_CT_QUAD n+1 ];

% Maximize problem and return the objective value
[rcode,res] = mosekopt('maximize echo(0)', prob, []);
x = res.sol.itr.xx;    
er = mu'*x;

portfolio_2_frontier.m

Listing 17.20 portfolio_2_frontier.m Click here to download.
%%
%  File : portfolio_2_frontier.m
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  Description :  Implements a basic portfolio optimization model.
%                 Determines points on the efficient frontier.
% 
%%

function portfolio_2_frontier()

n      = 3;
w      = 1.0;   
mu     = [0.1073,0.0737,0.0627]';
x0     = [0.0,0.0,0.0]';
GT     = [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238  ; ...
           0.0              , 0.102863378954911  , -0.00222873156550421 ; ...
           0.0              , 0.0                ,  0.0338148677744977 ];

% Some predefined alphas are chosen
alphas = 0.6:0.1:15.0;
front = EfficientFrontier(n,mu,GT,x0,w,alphas);
disp(sprintf('\nEfficient frontier'));
disp(sprintf('%-12s  %-12s  %-12s', 'alpha', 'return', 'risk'));
for p = front.'
    disp(sprintf('%-12.4f  %-12.4e  %-12.4e', p(1), p(2), p(3)));
end
plot(front(:,3),front(:,2));


%
%    Purpose:
%        Computes several portfolios on the optimal portfolios by 
%
%            for alpha in alphas: 
%                maximize   expected return - alpha * variance
%                subject to the constraints  
%        
%    Input:
%        n: Number of assets
%        mu: An n dimmensional vector of expected returns
%        GT: A matrix with n columns so (GT')*GT  = covariance matrix
%        x0: Initial holdings 
%        w: Initial cash holding
%        alphas: List of the alphas
%                    
%    Output:
%        The efficient frontier as list of tuples (alpha, expected return, variance)
%
function frontier = EfficientFrontier(n,mu,GT,x0,w,alphas)

frontier = [];
[rcode, res] = mosekopt('symbcon');
prob = [];

% The budget constraint in terms of variables [x; s]
prob.a = [ones(1,n), 0.0];
prob.blc = w + sum(x0);
prob.buc = w + sum(x0);

% No shortselling
prob.blx = [zeros(n,1); -inf];
prob.bux = inf*ones(n+1,1); 

% An affine conic constraint: [s, 0.5, GT*x] in rotated quadratic cone
% In matrix form
% [ 0  1]  [ x ]     [ 0   ]
% [ 0  0]  [   ]  +  [ 0.5 ]  \in Q_r
% [ GT 0]  [ s ]     [ 0   ]  
prob.f = sparse([ [zeros(1,n), 1.0]; zeros(1, n+1); [GT, zeros(n,1)] ]);
prob.g = [ 0; 0.5; zeros(n, 1) ]
prob.cones = [ res.symbcon.MSK_CT_RQUAD n+2 ];

for alpha = alphas
    % Objective mu'*x - alpha*s
    prob.c = [mu; -alpha];

    [rcode,res] = mosekopt('maximize echo(0)',prob,[]);
    x = res.sol.itr.xx(1:n);
    s = res.sol.itr.xx(n+1);    
       
    frontier = [frontier; [alpha, mu'*x, s] ];
end

portfolio_3_impact.m

Listing 17.21 portfolio_3_impact.m Click here to download.
%%
%  File : portfolio_3_impact.m
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  Description :  Implements a basic portfolio optimization model
%                 with transaction costs of order x^(3/2).
% 
%%

function portfolio_3_impact()

n      = 3;
w      = 1.0;   
mu     = [0.1073,0.0737,0.0627]';
x0     = [0.0,0.0,0.0]';
gammas = [0.035,0.040,0.050,0.060,0.070,0.080 ,0.090]';
GT     = [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238  ; ...
           0.0              , 0.102863378954911  , -0.00222873156550421 ; ...
           0.0              , 0.0                ,  0.0338148677744977 ];

% Somewhat arbirtrary choice of m
m = 1.0e-2*ones(n,1);
gamma = gammas(1);
[x, t] = MarkowitzWithMarketImpact(n,mu,GT,x0,w,gamma,m);
disp(sprintf('\nMarkowitz portfolio optimization with market impact cost'));
disp(sprintf('Expected return: %.4e Std. deviation: %.4e Market impact cost: %.4e', mu'*x, gamma, m'*t));

%
%        Description:
%            Extends the basic Markowitz model with a market cost term.
%
%        Input:
%            n: Number of assets
%            mu: An n dimmensional vector of expected returns
%            GT: A matrix with n columns so (GT')*GT  = covariance matrix
%            x0: Initial holdings 
%            w: Initial cash holding
%            gamma: Maximum risk (=std. dev) accepted
%            m: It is assumed that  market impact cost for the j'th asset is
%               m_j|x_j-x0_j|^3/2
%
%        Output:
%           Optimal expected return and the optimal portfolio     
%
function [x, t] = MarkowitzWithMarketImpact(n,mu,GT,x0,w,gamma,m)

[rcode, res] = mosekopt('symbcon');

% unrolled variable ordered as (x, t)
prob = [];
prob.c = [mu; zeros(n,1)];

In = speye(n);
On = sparse([],[],[],n,n);

% Linear part
% [ e' m' ]  * [ x; t ]  =   w + e'*x0
prob.a   = [ ones(1,n), m' ]; 
prob.blc = [ w + sum(x0) ];
prob.buc = [ w + sum(x0) ];

% No shortselling and no other bounds
prob.blx = [ zeros(n,1); -inf*ones(n,1) ];
prob.bux = inf*ones(2*n,1);

% Affine conic constraints representing [ gamma, GT*x ] in quadratic cone
prob.f = sparse([ zeros(1,2*n); [GT On] ]);
prob.g = [gamma; zeros(n,1)];
prob.cones = [ res.symbcon.MSK_CT_QUAD n+1 ];

% Extend the affine conic constraints
% with power cones representing t(i) >= |x(i)-x0(i)|^1.5
fi = [];
fj = [];
g  = [];
fv = repmat([1; 1], n, 1);
for k=1:n
    fi = [fi; 3*k-2; 3*k];
    fj = [fj; n+k; k];
    g  = [g; 0; 1; -x0(k)];
end
prob.f = [prob.f; sparse(fi, fj, fv)];
prob.g = [prob.g; g];
prob.cones = [prob.cones repmat([res.symbcon.MSK_CT_PPOW, 3, 2, 2.0, 1.0], 1, n) ];

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

x = res.sol.itr.xx(1:n);
t = res.sol.itr.xx(n+(1:n));

portfolio_4_transcost.m

Listing 17.22 portfolio_4_transcost.m Click here to download.
%%
%  File : portfolio_4_transcost.m
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  Description :  Implements a basic portfolio optimization model
%                 with fixed setup costs and transaction costs
%                 as a mixed-integer problem.
% 
%%

function portfolio_4_transcost()

n      = 3;
w      = 1.0;   
mu     = [0.1073,0.0737,0.0627]';
x0     = [0.0,0.0,0.0]';
gammas = [0.035,0.040,0.050,0.060,0.070,0.080 ,0.090]';
GT     = [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238  ; ...
           0.0              , 0.102863378954911  , -0.00222873156550421 ; ...
           0.0              , 0.0                ,  0.0338148677744977 ];

f = 0.01*ones(n,1);
g = 0.001*ones(n,1);
gamma = gammas(1);
[x, z, y] = MarkowitzWithTransactionsCost(n,mu,GT,x0,w,gamma,f,g);
disp(sprintf('\nMarkowitz portfolio optimization with transactions cost'))
disp(sprintf('Expected return: %.4e Std. deviation: %.4e Transactions cost: %.4e', ...
                 mu'*x,gamma,f'*y+g'*z))

%
%        Description:
%            Extends the basic Markowitz model with a market cost term.
%
%        Input:
%            n: Number of assets
%            mu: An n dimmensional vector of expected returns
%            GT: A matrix with n columns so (GT')*GT  = covariance matrix
%            x0: Initial holdings 
%            w: Initial cash holding
%            gamma: Maximum risk (=std. dev) accepted
%            f: If asset j is traded then a fixed cost f_j must be paid
%            g: If asset j is traded then a cost g_j must be paid for each unit traded
%
%        Output:
%           Optimal expected return and the optimal portfolio     
%
function [x, z, y] = MarkowitzWithTransactionsCost(n,mu,GT,x0,w,gamma,f,g)

[rcode, res] = mosekopt('symbcon');

% Upper bound on the traded amount
u = w+sum(x0);

% unrolled variable ordered as (x, z, y)
prob = [];
prob.c = [mu; zeros(2*n,1)];
In = speye(n);
On = sparse([],[],[],n,n);

% Linear constraints
% [ e'  g'  f' ]   [ x ]  =   w + e'*x0
% [ I  -I   0  ] * [ z ]  <=  x0
% [ I   I   0  ]   [ y ]  >=  x0
% [ 0   I  -U  ]          <=  0
prob.a   = [ [ones(1,n), g', f']; In -In On; In In On; On In -u*In ]; 
prob.blc = [ w + sum(x0); -inf*ones(n,1); x0; -inf*ones(n,1) ];
prob.buc = [ w + sum(x0); x0; inf*ones(n,1); zeros(n,1) ];

% No shortselling and the linear bound 0 <= y <= 1 
prob.blx = [ zeros(n,1); -inf*ones(n,1); zeros(n,1) ];
prob.bux = [ inf*ones(2*n,1); ones(n,1) ];

% Affine conic constraints representing [ gamma, GT*x ] in quadratic cone
prob.f = sparse([ zeros(1,3*n); [GT On On];  ]);
prob.g = [gamma; zeros(n,1)];
prob.cones = [ res.symbcon.MSK_CT_QUAD n+1 ];

% Demand y to be integer (hence binary)
prob.ints.sub = 2*n+(1:n);

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

x = res.sol.int.xx(1:n);
z = res.sol.int.xx(n+(1:n));
y = res.sol.int.xx(2*n+(1:n));

portfolio_5_card.m

Listing 17.23 portfolio_5_card.m Click here to download.
%%
%  File : portfolio_5_card.m
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  Description :  Implements a basic portfolio optimization model
%                 with cardinality constraints on number of assets traded.
% 
%%
function portfolio_5_card()

n      = 3;
w      = 1.0;   
mu     = [0.1073,0.0737,0.0627]';
x0     = [0.0,0.0,0.0]';
gammas = [0.035,0.040,0.050,0.060,0.070,0.080 ,0.090]';
GT     = [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238  ; ...
           0.0              , 0.102863378954911  , -0.00222873156550421 ; ...
           0.0              , 0.0                ,  0.0338148677744977 ];

disp(sprintf('\nMarkowitz portfolio optimization with cardinality constraints'));
for k=1:n
    x = MarkowitzWithCardinality(n,mu,GT,x0,w,gammas(1),k);
    disp(sprintf('Bound: %d   Expected return: %.4e  Solution:', k, mu'*x));
    disp(x);
end

%
%    Description:
%        Extends the basic Markowitz model with cardinality constraints.
%
%    Input:
%        n: Number of assets
%        mu: An n dimensional vector of expected returns
%        GT: A matrix with n columns so (GT')*GT  = covariance matrix
%        x0: Initial holdings 
%        w: Initial cash holding
%        gamma: Maximum risk (=std. dev) accepted
%        k: Maximum number of assets on which we allow to change position.
%
%    Output:
%       Optimal expected return and the optimal portfolio     
%
function x = MarkowitzWithCardinality(n,mu,GT,x0,w,gamma,k)

[rcode, res] = mosekopt('symbcon');

% Upper bound on the traded amount
u = w+sum(x0);

% unrolled variable ordered as (x, z, y)
prob = [];
prob.c = [mu; zeros(2*n,1)];
In = speye(n);
On = sparse([],[],[],n,n);

% Linear constraints
% [ e'  0   0  ]           =   w + e'*x0
% [ I  -I   0  ]   [ x ]  <=  x0
% [ I   I   0  ] * [ z ]  >=  x0
% [ 0   I  -U  ]   [ y ]  <=  0
% [ 0   0   e' ]          <=  k
prob.a   = [ [ones(1,n), zeros(1,2*n)]; In -In On; In In On; On In -u*In; zeros(1,2*n) ones(1,n) ]; 
prob.blc = [ w + sum(x0); -inf*ones(n,1); x0; -inf*ones(n,1); 0 ];
prob.buc = [ w + sum(x0); x0; inf*ones(n,1); zeros(n,1); k ];

% No shortselling and the linear bound 0 <= y <= 1 
prob.blx = [ zeros(n,1); -inf*ones(n,1); zeros(n,1) ];
prob.bux = [ inf*ones(2*n,1); ones(n,1) ];

% Affine conic constraints representing [ gamma, GT*x ] in quadratic cone
prob.f = sparse([ zeros(1,3*n); [GT On On];  ]);
prob.g = [gamma; zeros(n,1)];
prob.cones = [ res.symbcon.MSK_CT_QUAD n+1 ];

% Demand y to be integer (hence binary)
prob.ints.sub = 2*n+(1:n);

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

x = res.sol.int.xx(1:n);

pow1.m

Listing 17.24 pow1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      pow1.m
%
%  Purpose:   To demonstrate how to solve the problem
%
%     maximize x^0.2*y^0.8 + z^0.4 - x
%           st x + y + 0.5z = 2
%              x,y,z >= 0

function pow1()

clear prob;

[r, res] = mosekopt('symbcon');
% Specify the non-conic part of the problem.

prob.c   = [-1 0 0 1 1 0];
prob.a   = [1 1 0.5 0 0 0];
prob.blc = [2.0];
prob.buc = [2.0];
prob.blx = [-inf -inf -inf -inf -inf 1.0];
prob.bux = [ inf  inf  inf  inf  inf 1.0];

% Specify the cones.
prob.cones.type   = [res.symbcon.MSK_CT_PPOW res.symbcon.MSK_CT_PPOW];
prob.cones.conepar= [0.2 0.4];
prob.cones.sub    = [1 2 4 3 6 5];
prob.cones.subptr = [1 4];
% The field 'type' specifies the cone types, in this case power cones.
%
% The fields 'sub' and 'subptr' specify the members of the cones,
% i.e., the above definitions imply that 
%   (x(1), x(2), x(4))  and (x(3), x(6), x(5)) 
% are cones.
%
% The field 'conepar' specifies the alpha cone parameters (exponents)

% Optimize the problem. 

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

% Display the primal solution.

res.sol.itr.xx'

qcqo1.m

Listing 17.25 qcqo1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      qcqo1.m
%
%  Purpose : Demonstrates a simple quadratically constrained quadratic problem 
%            minimize  x_1^2 + 0.1 x_2^2 +  x_3^2 - x_1 x_3 - x_2
%            s.t 1 <=  x_1 + x_2 + x_3 - x_1^2 - x_2^2 - 0.1 x_3^2 + 0.2 x_1 x_3
%            x >= 0
%
function qcqo1()
clear prob;

% Specify the linear objective terms.
prob.c      = [0, -1, 0];

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

% Specify the quadratic terms of the objective.
prob.qosubi = [1     2    3    3  ]';
prob.qosubj = [1     2    3    1  ]';
prob.qoval  = [2.0   0.2  2.0 -1.0]';

% Specify the linear constraint matrix
prob.a      = [1 1 1];

% Specify the lower bounds
prob.blc    = [1];
prob.blx    = zeros(3,1);

[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));

qo1.m

Listing 17.26 qo1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      qo1.m
%
%  Purpose : Demonstrates a simple quadratic problem 
%
function qo1()

% 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

qo2.m

Listing 17.27 qo2.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      qo2.m
%
%  Purpose : Demonstrates a simple quadratic problem 
%
function qo2()

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'

reoptimization.m

Listing 17.28 reoptimization.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      reoptimization.m
%
%  Purpose : Demonstrates how to modify and re-optimize a linear problem
%

function reoptimization()
clear prob;

% Specify the c vector.
prob.c  = [1.5 2.5 3.0]';

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

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

% Specify lower bounds of the constraints.
prob.blc = [-inf -inf -inf]';

% Specify  upper bounds of the constraints.
prob.buc = [100000 50000 60000]';

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

% Specify upper bounds of the variables.
prob.bux = [inf inf inf]';

% Perform the optimization.
param.MSK_IPAR_OPTIMIZER = 'MSK_OPTIMIZER_FREE_SIMPLEX'; 
[r,res] = mosekopt('maximize',prob,param); 

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change a coefficient
prob.a(1,1) = 3.0;

% Reoptimize with changed coefficient
% Use previous solution to perform very simple hot-start.
% This part can be skipped, but then the optimizer will start
% from scratch on the new problem, i.e. without any hot-start.
prob.sol = [];
prob.sol.bas = res.sol.bas;			
[r,res] = mosekopt('maximize',prob,param);
res.sol.bas.xx

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add a variable
prob.c       = [prob.c;1.0];
prob.a       = [prob.a,sparse([4.0 0.0 1.0]')];
prob.blx     = [prob.blx; 0.0];
prob.bux     = [prob.bux; inf];

% Reoptimize with a new variable and hot-start
% All parts of the solution must be extended to the new dimensions.
prob.sol = [];
prob.sol.bas = res.sol.bas;			
prob.sol.bas.xx  = [prob.sol.bas.xx; 0.0];
prob.sol.bas.slx = [prob.sol.bas.slx; 0.0];
prob.sol.bas.sux = [prob.sol.bas.sux; 0.0];
prob.sol.bas.skx = [prob.sol.bas.skx; 'UN'];
[r,res] = mosekopt('maximize',prob,param);
res.sol.bas.xx

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add a constraint
prob.a       = [prob.a;sparse([1.0 2.0 1.0 1.0])]; 
prob.blc     = [prob.blc; -inf]; 
prob.buc     = [prob.buc; 30000.0]; 

% Reoptimize with a new variable and hot-start
prob.sol = [];
prob.sol.bas = res.sol.bas;			
prob.sol.bas.y   = [prob.sol.bas.y; 0.0];
prob.sol.bas.xc  = [prob.sol.bas.xc; 0.0];
prob.sol.bas.slc = [prob.sol.bas.slc; 0.0];
prob.sol.bas.suc = [prob.sol.bas.suc; 0.0];
prob.sol.bas.skc = [prob.sol.bas.skc; 'UN'];
[r,res] = mosekopt('maximize',prob,param);
res.sol.bas.xx


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change constraint bounds
prob.buc     = [80000 40000 50000 22000]';
prob.sol     = res.sol;
[r,res] = mosekopt('maximize',prob,param);
res.sol.bas.xx

response.m

Listing 17.29 response.m Click here to download.
%%
%  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File:      response.m
%
%  Purpose:   This example demonstrates proper response handling
%             for problems solved with the interior-point optimizers.
%

function response(inputfile)

  cmd = sprintf('read(%s)', inputfile)
  % In this example we read the problem from file
  [r, res] = mosekopt(cmd)

  % Read was successful
  if strcmp(res.rcodestr, 'MSK_RES_OK')

      prob = res.prob;
      param = []

      % (Optionally) Uncomment the next line to get solution status Unknown
      % param.MSK_IPAR_INTPNT_MAX_ITERATIONS = 1

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

      % Expected result: The solution status of the interior-point solution is optimal.

      % Check if we have non-error termination code or OK
      if isempty(strfind(res.rcodestr, 'MSK_RES_ERR'))
      
          solsta = strcat('MSK_SOL_STA_', res.sol.itr.solsta)

          if strcmp(solsta, 'MSK_SOL_STA_OPTIMAL')
              fprintf('An optimal interior-point solution is located:\n');
              res.sol.itr.xx
              
          elseif strcmp(solsta, 'MSK_SOL_STA_DUAL_INFEASIBLE_CER')
              fprintf('Dual infeasibility certificate found.');

          elseif strcmp(solsta, 'MSK_SOL_STA_PRIMAL_INFEASIBLE_CER')
              fprintf('Primal infeasibility certificate found.');

          elseif strcmp(solsta, 'MSK_SOL_STA_UNKNOWN') 
              % The solutions status is unknown. The termination code 
              % indicates why the optimizer terminated prematurely. 
              fprintf('The solution status is unknown.\n');
              fprintf('Termination code: %s (%d) %s.\n', res.rcodestr, res.rcode, res.rmsg);  
          else
            fprintf('An unexpected solution status is obtained.');
          end
      
      else
          fprintf('Error during optimization: %s (%d) %s.\n', res.rcodestr, res.rcode, res.rmsg);  
      end

  else
      fprintf('Could not read input file, error: %s (%d) %s.\n', res.rcodestr, res.rcode, res.rmsg)
  end

end

rlo1.m

Listing 17.30 rlo1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      rlo1.m
%
%  Purpose : First part of a robust optimization example (see rlo2.m)
%
function rlo1()

prob.c   = [-100;-199.9;6200-700;6900-800];
prob.a   = sparse([0.01,0.02,-0.500,-0.600;1,1,0,0;
                   0,0,90.0,100.0;0,0,40.0,50.0;100.0,199.9,700,800]);
prob.blc = [0;-inf;-inf;-inf;-inf];
prob.buc = [inf;1000;2000;800;100000];
prob.blx = [0;0;0;0];
prob.bux = [inf;inf;inf;inf];
[r,res]  = mosekopt('maximize',prob);
xx       = res.sol.itr.xx;
RawI     = xx(1);
RawII    = xx(2);
DrugI    = xx(3);
DrugII   = xx(4);

disp(sprintf('*** Optimal value: %8.3f',prob.c'*xx));
disp('*** Optimal solution:');
disp(sprintf('RawI:   %8.3f',RawI));
disp(sprintf('RawII:  %8.3f',RawII));
disp(sprintf('DrugI:  %8.3f',DrugI));
disp(sprintf('DrugII: %8.3f',DrugII));

rlo2.m

Listing 17.31 rlo2.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      rlo2.m
%
%  Purpose :   Solves the problem:
%
%              Maximize t 
%              subject to
%                        t <= sum(delta(j)*x(j)) -Omega*z,
%                        y(j) = sigma(j)*x(j), j=1,...,n,
%                        sum(x(j)) = 1,
%                        || y || <= z,
%

function rlo2(n, Omega, draw)

n = str2num(n)
Omega = str2num(Omega)
draw

[r, res] = mosekopt('symbcon echo(0)');
sym = res.symbcon;

% Set nominal returns and volatilities
delta = (0.96/(n-1))*[0:1:n-1]+1.04;
sigma = (1.152/(n-1))*[0:1:n-1];

% Set mosekopt description of the problem
prob.c             = [1;zeros(n+1,1)];
A                  = [-1, delta,     -Omega; ...
                       0, ones(1,n), 0];
prob.a             = sparse(A);
prob.blc           = [0;1];
prob.buc           = [inf;1];
prob.blx           = [-inf;zeros(n,1);0];
prob.bux           = inf*ones(n+2,1);
F                  = [zeros(1,n+1), 1; ...
                      zeros(n,1), diag(sigma), zeros(n,1)];
prob.f             = sparse(F);
prob.cones         = [ sym.MSK_CT_QUAD n+1 ];

% Run mosekopt
[r,res]=mosekopt('maximize echo(3)',prob);

xx = res.sol.itr.xx;
t  = xx(1);

disp(sprintf('Robust optimal value: %5.4f',t));
x = max(xx(2:1+n),zeros(n,1));

if draw == true 
    % Display the solution
    plot([1:1:n],x,'-m');
    grid on;
    
    disp('Press <Enter> to run simulations');
    pause
    
    % Run simulations
    
    Nsim = 10000;
    out  = zeros(Nsim,1);
    for i=1:Nsim,
        returns  = delta+(2*rand(1,n)-1).*sigma;
        out(i)   = returns*x;
    end;
    disp(sprintf('Actual returns over %d simulations:',Nsim));
    disp(sprintf('Min=%5.4f Mean=%5.4f Max=%5.4f StD=%5.2f',...
                 min(out),mean(out),max(out),std(out)));
    hist(out);
end

sdo1.m

Listing 17.32 sdo1.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      sdo1.m
%
%  Purpose :   Solves the mixed semidefinite and conic quadratic optimization problem
%
%                 minimize    Tr [2, 1, 0; 1, 2, 1; 0, 1, 2]*X + x(1)
%
%                 subject to  Tr [1, 0, 0; 0, 1, 0; 0, 0, 1]*X + x(1)               = 1
%                             Tr [1, 1, 1; 1, 1, 1; 1, 1, 1]*X        + x(2) + x(3) = 0.5
%                             X>=0,  x(1) >= sqrt(x(2)^2 + x(3)^2)
%


function sdo1()
[r, res] = mosekopt('symbcon');

prob.c         = [1, 0, 0];

prob.bardim    = [3];
prob.barc.subj = [1, 1, 1, 1, 1];
prob.barc.subk = [1, 2, 2, 3, 3];
prob.barc.subl = [1, 1, 2, 2, 3];
prob.barc.val  = [2.0, 1.0, 2.0, 1.0, 2.0];

prob.blc = [1, 0.5];
prob.buc = [1, 0.5];

% It is a good practice to provide the correct 
% dimmension of A as the last two arguments
% because it facilitates better error checking.
prob.a         = sparse([1, 2, 2], [1, 2, 3], [1, 1, 1], 2, 3);
prob.bara.subi = [1, 1, 1, 2, 2, 2, 2, 2, 2];
prob.bara.subj = [1, 1, 1, 1, 1, 1, 1, 1, 1];
prob.bara.subk = [1, 2, 3, 1, 2, 3, 2, 3, 3];
prob.bara.subl = [1, 2, 3, 1, 1, 1, 2, 2, 3];
prob.bara.val  = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];

prob.cones.type   = [res.symbcon.MSK_CT_QUAD];
prob.cones.sub    = [1, 2, 3];
prob.cones.subptr = [1];

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

X = zeros(3);
X([1,2,3,5,6,9]) = res.sol.itr.barx;
X = X + tril(X,-1)';

x = res.sol.itr.xx;

sdo2.m

Listing 17.33 sdo2.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      sdo2.m
%
%  Purpose :   Solves the semidefinite problem with two symmetric variables:
%
%                 min   <C1,X1> + <C2,X2>
%                 st.   <A1,X1> + <A2,X2> = b
%                             (X2)_{1,2} <= k
%                
%                 where X1, X2 are symmetric positive semidefinite,
%
%                 C1, C2, A1, A2 are assumed to be constant symmetric matrices,
%                 and b, k are constants.
                
function sdo2()

% Sample data
C1 = [ 1 0 0; 0 0 0; 0 0 6 ];
A1 = [ 1 0 1; 0 0 0; 1 0 2 ];
C2 = [ 1 -3 0 0; -3 2 0 0; 0 0 1 0; 0 0 0 0 ];
A2 = [ 0 1 0 0; 1 -1 0 0; 0 0 0 0; 0 0 0 -3 ];
b = 23;
k = -3;

% The scalar part, as in linear optimization examples
prob.c = [];
prob.a = sparse([], [], [], 2, 0);          % 2 constraints, no scalar variables
prob.blc = [b -inf];                        % Bounds
prob.buc = [b k];

% Dimensions of PSD variables
prob.bardim = [3, 4];

% Coefficients in the objective
[r1,c1,v1] = find(tril(C1));
[r2,c2,v2] = find(tril(C2));

prob.barc.subj = [repmat(1,length(v1),1);                % Which PSD variable (j)
                  repmat(2,length(v2),1)];           
prob.barc.subk = [r1; r2];                               % Which matrix entry and value ((k,l)->v)
prob.barc.subl = [c1; c2];
prob.barc.val = [v1; v2];

% Coefficients in the constraints
[r1,c1,v1] = find(tril(A1));
[r2,c2,v2] = find(tril(A2));

prob.bara.subi = [ones(length(v1)+length(v2),1);         % Which constraint (i)
                  2];     
prob.bara.subj = [repmat(1,length(v1),1); 
                  repmat(2,length(v2),1); 
                  2];                                    % Which PSD variable (j)
prob.bara.subk = [r1; r2; 2];                            % Which matrix entry and value ((k,l)->v)
prob.bara.subl = [c1; c2; 1];
prob.bara.val = [v1; v2; 0.5];

% Solve with log output
[r, res] = mosekopt('write(test.ptf) minimize echo(10)', prob);

% Retrieve the result assuming primal and dual feasible
X1 = zeros(3);
X1([1,2,3,5,6,9]) = res.sol.itr.barx(1:6);
X1 = X1 + tril(X1,-1)';

X2 = zeros(4);
X2([1,2,3,4,6,7,8,11,12,16]) = res.sol.itr.barx(7:16);
X2 = X2 + tril(X2,-1)';

X1
X2

sensitivity.m

Listing 17.34 sensitivity.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      sensitivity.m
%
%  Purpose:   To demonstrate how to perform sensitivity
%             analysis from the MOSEK Toolbox on a small problem:
%
%  minimize
%
%  obj: +1 x11 + 2 x12 + 5 x23 + 2 x24 + 1 x31 + 2 x33 + 1 x34
%  st
%  c1:   +  x11 +   x12                                          <= 400
%  c2:                  +   x23 +   x24                          <= 1200
%  c3:                                  +   x31 +   x33 +   x34  <= 1000
%  c4:   +  x11                         +   x31                   = 800
%  c5:          +   x12                                           = 100
%  c6:                  +   x23                 +   x33           = 500
%  c7:                          +   x24                 +   x34   = 500
%
function sensitivity()

clear prob;

% Obtain all symbolic constants
% defined by MOSEK.
[r,res]  = mosekopt('symbcon');
sc       = res.symbcon;

prob.blc = [-Inf, -Inf, -Inf, 800,100,500,500];
prob.buc = [ 400, 1200, 1000, 800,100,500,500];
prob.c   = [1.0,2.0,5.0,2.0,1.0,2.0,1.0]';
prob.blx = [0.0,0.0,0.0,0.0,0.0,0.0,0.0];
prob.bux = [Inf,Inf,Inf,Inf, Inf,Inf,Inf];

subi     = [  1,  1,  2,  2,  3,  3,  3,  4,  4,  5,  6,  6,  7,  7];
subj     = [  1,  2,  3,  4,  5,  6,  7,  1,  5,  6,  3,  6,  4,  7];
val      = [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0];

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

% analyse upper bound 1:7
prob.prisen.cons.subl = [];
prob.prisen.cons.subu = [1:7];
% analyse lower bound on variables 1:7
prob.prisen.vars.subl = [1:7];
prob.prisen.vars.subu = [];
% analyse coeficient 1:7
prob.duasen.sub = [1:7];
[r,res] = mosekopt('minimize echo(0)',prob);

%Print results

fprintf('\nBasis sensitivity results:\n')
fprintf('\nSensitivity for bounds on constraints:\n')
for i = 1:length(prob.prisen.cons.subl)
  fprintf (...
  'con = %d, beta_1 = %.1f, beta_2 = %.1f, delta_1 = %.1f,delta_2 = %.1f\n', ...
  prob.prisen.cons.subl(i),res.prisen.cons.lr_bl(i), ...
  res.prisen.cons.rr_bl(i),...
  res.prisen.cons.ls_bl(i),...
  res.prisen.cons.rs_bl(i));
end

for i = 1:length(prob.prisen.cons.subu)
  fprintf (...
  'con = %d, beta_1 = %.1f, beta_2 = %.1f, delta_1 = %.1f,delta_2 = %.1f\n', ...
  prob.prisen.cons.subu(i),res.prisen.cons.lr_bu(i), ...
  res.prisen.cons.rr_bu(i),...
  res.prisen.cons.ls_bu(i),...
  res.prisen.cons.rs_bu(i));
end
fprintf('Sensitivity for bounds on variables:\n')
for i = 1:length(prob.prisen.vars.subl)
fprintf (...
'var = %d, beta_1 = %.1f, beta_2 = %.1f, delta_1 = %.1f,delta_2 = %.1f\n', ...
 prob.prisen.vars.subl(i),res.prisen.vars.lr_bl(i), ...
 res.prisen.vars.rr_bl(i),...
 res.prisen.vars.ls_bl(i),...
 res.prisen.vars.rs_bl(i));
end

for i = 1:length(prob.prisen.vars.subu)
  fprintf (...
  'var = %d, beta_1 = %.1f, beta_2 = %.1f, delta_1 = %.1f,delta_2 = %.1f\n', ...
  prob.prisen.vars.subu(i),res.prisen.vars.lr_bu(i), ...
  res.prisen.vars.rr_bu(i),...
  res.prisen.vars.ls_bu(i),...
  res.prisen.vars.rs_bu(i));
end

fprintf('Sensitivity for coefficients in objective:\n')
for i = 1:length(prob.duasen.sub)
  fprintf (...
  'var = %d, beta_1 = %.1f, beta_2 = %.1f, delta_1 = %.1f,delta_2 = %.1f\n', ...
  prob.duasen.sub(i),res.duasen.lr_c(i), ...
  res.duasen.rr_c(i),...
  res.duasen.ls_c(i),...
  res.duasen.rs_c(i));
end

sensitivity2.m

Listing 17.35 sensitivity2.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      sensitivity2.m
%

function sensitivity2()
% Setup problem data.

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  );

simple.m

Listing 17.36 simple.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      simple.m
%
%  Purpose :   To demonstrate how solve a problem
%              read from file.
%

function simple(inputfile, solfile)

cmd      = sprintf('read(%s)', inputfile)
% Read the problem from file
[rcode, res] = mosekopt(cmd)

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

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



end

solutionquality.m

Listing 17.37 solutionquality.m Click here to download.
%
%  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
%  File :      solutionquality.m
%
%  Purpose :   To demonstrate how to examine the quality of a solution. 
%
 
function solutionquality(data)

    cmd      = sprintf('read(%s)',data)
    % Read the problem from file
    [r, res] = mosekopt(cmd)

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

    solsta = strcat('MSK_SOL_STA_', res.sol.itr.solsta);
    
    if strcmp(solsta, 'MSK_SOL_STA_OPTIMAL')

      sol = res.sol.itr
      primalobj= sol.pobjval
      dualobj= sol.dobjval
      
      abs_obj_gap     = abs(dualobj - primalobj);
      rel_obj_gap     = abs_obj_gap/(1.0 + min( abs(primalobj), abs(dualobj)));

      
      % Assume the application needs the solution to be within
      % 1e-6 optimality in an absolute sense. Another approach
      % would be looking at the relative objective gap */
      
      fprintf('\n\n');
      fprintf('Customized solution information.\n');
      fprintf('  Absolute objective gap: %e\n',abs_obj_gap);
      fprintf('  Relative objective gap: %e\n',rel_obj_gap);
          
      if ( rel_obj_gap>1e-6 )
          fprintf('Warning: The relative objective gap is LARGE.\n');
          % Perhaps we reject the solution here
      else            
          % We accept solution
          res.sol.itr.xx
      end
          
   
    elseif strcmp(solsta, 'MSK_SOL_STA_DUAL_INFEASIBLE_CER') || ...
           strcmp(solsta, 'MSK_SOL_STA_PRIMAL_INFEASIBLE_CER') 
        fprintf('Primal or dual infeasibility certificate found.\n');
   
    elseif strcmp(solsta,  'MSK_SOL_STA_UNKNOWN')
        fprintf('The status of the solution is unknown.\n');

    else
        fprintf('Other solution status');
    end

end