17 List of examples¶
List of examples shipped in the distribution of Optimization Toolbox for MATLAB:
File |
Description |
|---|---|
Advanced simplex hot-start examples |
|
A simple problem using affine conic constraints |
|
A simple problem using affine conic constraints |
|
An example of data/progress callback |
|
Log handler definition for |
|
A simple conic exponential problem |
|
A simple conic quadratic problem |
|
A simple example of how to repair an infeasible problem |
|
A simple geometric program (GP) in conic form |
|
A Hello World example |
|
A simple linear problem using msklpopt |
|
A simple linear problem using mosekopt |
|
A simple linear problem using linprog |
|
A simple mixed-integer conic problem |
|
A simple mixed-integer linear problem |
|
Smallest disk covering a subset of points (MICQO) |
|
A simple mixed-integer linear problem with an initial guess |
|
Demonstrates least squares and other norm minimization problems |
|
Uses MOSEK OptServer to solve an optimization problem synchronously |
|
Shows how to set optimizer parameters and read information items |
|
Shows how to obtain and analyze a primal infeasibility certificate |
|
Portfolio optimization - basic Markowitz model |
|
Portfolio optimization - efficient frontier |
|
Portfolio optimization - market impact costs |
|
Portfolio optimization - transaction costs |
|
Portfolio optimization - cardinality constraints |
|
Portfolio optimization - factor model |
|
A simple power cone problem |
|
A simple quadratically constrained quadratic problem |
|
A simple quadratic problem |
|
A simple quadratic problem |
|
Demonstrate how to modify and re-optimize a linear problem |
|
Demonstrates proper response handling |
|
Robust linear optimization example, part 1 |
|
Robust linear optimization example, part 2 |
|
A simple semidefinite problem with one matrix variable and a quadratic cone |
|
A simple semidefinite problem with two matrix variables |
|
A simple semidefinite problem with an LMI using the SVEC domain. |
|
Sensitivity analysis performed on a small linear problem |
|
Sensitivity analysis performed on a small linear problem |
|
A simple I/O example: read problem from a file, solve and write solutions |
|
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
%%
% 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
%%
% 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_DOMAIN_QUADRATIC_CONE 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_DOMAIN_PRIMAL_POWER_CONE 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_DOMAIN_PRIMAL_POWER_CONE 3 2 1.0 3.0];
% All cones
prob.f = [FQ; FP1; FP2];
prob.g = [gQ; gP1; gP2];
prob.accs = [cQ cP1 cP2];
[r, res] = mosekopt('maximize', prob);
res.sol.itr.pobjval
res.sol.itr.xx(1:2)
affco2.m
%%
% 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.accs = [res.symbcon.MSK_DOMAIN_QUADRATIC_CONE n+1 repmat([res.symbcon.MSK_DOMAIN_PRIMAL_EXP_CONE 3], 1, n)];
% Solve
[r, res] = mosekopt('minimize', prob);
t = res.sol.itr.xx(1)
callback.m
%
% 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
% 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
%
% 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 affine conic constraint with one exponential cone.
prob.accs = [res.symbcon.MSK_DOMAIN_PRIMAL_EXP_CONE 3];
prob.f = speye(3);
% prob.accs the domain types, in this case a single exponential cone
% The matrix f is the ientity, meaning that
%
% I * x \in EXP
%
% which is exactly
%
% 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
%
% 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 as affine conic constraints.
% Two conic constrains: one with QUAD, one with RQUAD, both of dimension 3
prob.accs = [res.symbcon.MSK_DOMAIN_QUADRATIC_CONE 3 res.symbcon.MSK_DOMAIN_RQUADRATIC_CONE 3];
% The matrix such that f * x = [x(4), x(1), x(2), x(5), x(6), x(3)]
prob.f = sparse( 1:6, [4, 1, 2, 5, 6, 3], ones(1, 6) );
% That implies:
% (x(4), x(1), x(2)) \in QUAD_3
% (x(5), x(6), x(3)) \in RQUAD_3
% Optimize the problem.
[r,res]=mosekopt('minimize',prob);
% Display the primal solution.
res.sol.itr.xx'
feasrepairex1.m
%%
% 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
%
% 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 affine 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.accs = [ res.symbcon.MSK_DOMAIN_PRIMAL_EXP_CONE, 3, res.symbcon.MSK_DOMAIN_PRIMAL_EXP_CONE, 3 ];
% Optimize and print results
[r,res]=mosekopt('maximize',prob);
exp(res.sol.itr.xx(1:3))
helloworld.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: helloworld.m
%
% The most basic example of how to get started with MOSEK.
prob.a = sparse(0,1) % 0 linear constraints, 1 variable
prob.c = [1.0]' % Only objective coefficient
prob.blx= [2.0]' % Lower bound(s) on variable(s)
prob.bux= [3.0]' % Upper bound(s) on variable(s)
% Optimize
[r, res] = mosekopt('minimize', prob);
% Print answer
res.sol.itr.xx
lo1.m
%
% 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
%
% 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
%
% 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
%
% 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.accs = [symbcon.MSK_DOMAIN_QUADRATIC_CONE 3 symbcon.MSK_DOMAIN_PRIMAL_EXP_CONE 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 = [100, 9, -1];
% Optimize the problem.
[r,res] = mosekopt('minimize info',prob);
% The integer solution (x,y)
res.sol.int.xx(2:3)
milo1.m
%
% 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
mindisk.m
%%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File : mindisk.m
%
% Purpose: A mixed-integer conic quadratic
% problem.
%
% Given points (x(i),y(i)) and a bound L it finds
% the smallest disk which covers some L of the points.
%%
function mindisk()
% Coordinates of the given points (example)
x = [5.108736e+02 4.848650e+02 7.847473e+01 4.703155e+02 5.461198e+02 5.690343e+02 8.041608e+02 7.777201e+02 4.740465e+02];
y = [-1.906827e+02 -6.262071e+02 -6.383325e+02 -7.925310e+02 -3.453408e+02 -7.685773e+02 -9.449286e+02 -1.694813e+02 -7.486016e+02];
% How many points at least should be covered
L = 5;
% Big-M for the interger formulation
M = 3000;
% Initialize symbolic constants
[r,res] = mosekopt('symbcon');
sym = res.symbcon;
% Variables [R x y u1 ... un]
n = length(x);
prob = [];
prob.c = [1; zeros(n+2,1)];
% Constraint sum(u) >= L
prob.a = sparse([0 0 0 ones(1,n)]);
prob.blc = [L];
prob.buc = [inf];
% ui are binary
prob.ints.sub = 4:n+3;
prob.blx = [-inf;-inf;-inf; zeros(n,1)];
prob.bux = [ inf; inf; inf; ones(n,1)];
% A sequence of constraints
% (R + M(1-ui), x-xi, y-yi) \in Q3
prob.f=[];
prob.g=[];
for i=1:n
prob.f=[prob.f; eye(3) sparse(1, i, -M, 3, n)];
prob.g=[prob.g; M; -x(i); -y(i)];
end
prob.accs = repmat([sym.MSK_DOMAIN_QUADRATIC_CONE 3], 1,n);
% Solve
[r,res] = mosekopt('minimize', prob);
sol = res.sol.int.xx;
R = sol(1);
X = sol(2);
Y = sol(3);
% Optionally plot the result
%scatter(x, y, '.');
%hold on;
%fplot(@(t) X + R*cos(t), @(t) Y + R*sin(t))
%axis equal;
end
mioinitsol.m
%%
% 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]';
% Request constructing the solution from integer variable values
param.MSK_IPAR_MIO_CONSTRUCT_SOL = 1;
[r,res] = mosekopt('maximize info',prob,param);
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
%%
%
% 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()
% Change this to display the plots
usegraphics=false;
% 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);
if usegraphics
figure(1);
scatter(x1, y1);
hold on
end
% Least squares regression
x = norm_lse(xp, y1, zeros(0,6), []);
if usegraphics
plot(0:0.01:1, smallVander((0:0.01:1)', 5)*x);
end
% With ridge regularization
x = norm_lse_reg(xp, y1, zeros(0,6), [], 0.1);
if usegraphics
plot(0:0.01:1, smallVander((0:0.01:1)', 5)*x);
end
% Quadratic ridge
x = norm_lse_reg_quad(xp, y1, zeros(0,6), [], 0.2);
if usegraphics
plot(0:0.01:1, smallVander((0:0.01:1)', 5)*x);
end
% 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];
parts = [];
labels = [];
if usegraphics
figure(2);
scatter(x12, y12);
hold on
end
for p=[1.1, 2.0, 3.0, 6.0]
x = norm_p_norm(xp2, y12, zeros(0,6), [], p);
if usegraphics
parts = [parts; plot(0:0.01:1, smallVander((0:0.01:1)', 5)*x);];
labels = [labels; sprintf('%.1f', p)];
end
end
if usegraphics
legend(parts, labels);
end
end
% Least squares regression
% minimize \|Fx-g\|_2
function x = norm_lse(F,g,A,b)
clear prob;
[r, res] = mosekopt('symbcon echo(0)');
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.accs = [ res.symbcon.MSK_DOMAIN_QUADRATIC_CONE 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 echo(0)');
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.accs = [ res.symbcon.MSK_DOMAIN_QUADRATIC_CONE k+1 res.symbcon.MSK_DOMAIN_QUADRATIC_CONE 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 echo(0)');
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.accs = [ res.symbcon.MSK_DOMAIN_RQUADRATIC_CONE k+2 res.symbcon.MSK_DOMAIN_RQUADRATIC_CONE 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 echo(0)');
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.accs = [ res.symbcon.MSK_DOMAIN_QUADRATIC_CONE 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 echo(0)');
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.accs = [ repmat([res.symbcon.MSK_DOMAIN_PRIMAL_POWER_CONE, 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
%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File : opt_server_sync.m
%
% Purpose : Demonstrates solving a problem remotely using the OptServer.
%
% url should be 'http://server:port' or 'https://server:port'
% cert is the path to a TLS certificate, if using https
%
function opt_server_sync(inputfile, url, cert)
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
optserver.host = url;
param.MSK_SPAR_REMOTE_TLS_CERT_PATH = cert;
% Perform the optimization with full log output.
[r,res] = mosekopt(sprintf('%s echo(10)', prob.objsense), prob, param, [], optserver);
% Use the optimal x solution.
xx = res.sol.bas.xx;
parameters.m
%
% 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
pinfeas.m
%%
% File : pinfeas.m
%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% Purpose: Demonstrates how to fetch a primal infeasibility certificate
% for a linear problem
%%
function pinfeas()
% In this example we set up a simple problem
[prob] = testProblem();
% Perform the optimization.
[rcode, res] = mosekopt('minimize', prob);
% Check problem status
if strcmp(res.sol.itr.prosta, 'PRIMAL_INFEASIBLE')
% Set the tolerance at which we consider a dual value as essential
eps = 1e-7;
disp("Variable bounds important for infeasibility: ");
analyzeCertificate(res.sol.itr.slx, res.sol.itr.sux, eps);
disp("Constraint bounds important for infeasibility: ")
analyzeCertificate(res.sol.itr.slc, res.sol.itr.suc, eps);
else
disp("The problem is not primal infeasible, no certificate to show");
end
% Set up a simple linear problem from the manual for test purposes
function [prob] = testProblem()
prob = [];
prob.c = [1, 2, 5, 2, 1, 2, 1];
prob.a = sparse([1,1,2,2,3,3,3,4,4,5,6,6,7,7],...
[1,2,3,4,5,6,7,1,5,2,3,6,4,7],...
[1,1,1,1,1,1,1,1,1,1,1,1,1,1],...
7, 7);
prob.blc = [-inf, -inf, -inf, 1100, 200, 500, 500];
prob.buc = [200, 1000, 1000, 1100, 200, 500, 500];
prob.blx = [0, 0, 0, 0, 0, 0, 0];
prob.bux = repmat(inf, 1, 7);
prob
% Analyzes and prints infeasibility contributing elements
% sl - dual values for lower bounds
% su - dual values for upper bounds
% eps - tolerance for when a nunzero dual value is significant
function analyzeCertificate(sl, su, eps)
n = size(sl);
for i=1:n
if abs(sl(i)) > eps
disp(sprintf("#%d: lower, dual = %e", i, sl(i)));
end
if abs(su(i)) > eps
disp(sprintf("#%d: upper, dual = %e", i, su(i)));
end
end
portfolio_1_basic.m
%%
% 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 = 8;
w = 59.0;
mu = [0.07197349 0.15518171 0.17535435 0.0898094 0.42895777 0.39291844 0.32170722 0.18378628]';
x0 = [8.0 5.0 3.0 5.0 2.0 9.0 3.0 6.0]';
gamma = 36.0;
GT = [ 0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638; ...
0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506; ...
0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914; ...
0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742; ...
0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ; ...
0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187; ...
0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327; ...
0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ];
disp('Basic Markowitz portfolio optimization')
[er, x] = BasicMarkowitz(n,mu,GT,x0,w,gamma);
disp(sprintf('Expected return: %.4e Std. deviation: %.4e', er, gamma));
%
% 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, x] = BasicMarkowitz(n,mu,GT,x0,w,gamma)
[rcode, res] = mosekopt('symbcon echo(0)');
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.accs = [ res.symbcon.MSK_DOMAIN_QUADRATIC_CONE n+1 ];
% Maximize problem and return the objective value
[rcode,res] = mosekopt('maximize echo(0)', prob, []);
% Check if the interior point solution is an optimal point
% See https://docs.mosek.com/latest/toolbox/accessing-solution.html about handling solution statuses.
if startsWith(res.rcodestr, 'MSK_RES_ERR') | ~strcmp(res.sol.itr.solsta, 'OPTIMAL')
throw(MException('mosek:solsta', 'Not an optimal solution'));
end
x = res.sol.itr.xx;
er = mu'*x;
portfolio_2_frontier.m
%%
% 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()
% Change this to produce the plots
usegraphics = false;
n = 8;
w = 1.0;
mu = [0.07197349 0.15518171 0.17535435 0.0898094 0.42895777 0.39291844 0.32170722 0.18378628]';
x0 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]';
gamma = 36.0;
GT = [ 0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638; ...
0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506; ...
0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914; ...
0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742; ...
0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ; ...
0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187; ...
0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327; ...
0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ];
% Some predefined alphas are chosen
alphas = [0.0 0.01 0.1 0.25 0.30 0.35 0.4 0.45 0.5 0.75 1.0 1.5 2.0 3.0 10.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
if usegraphics
plot(front(:,3),front(:,2));
end
%
% 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 echo(0)');
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.accs = [ res.symbcon.MSK_DOMAIN_RQUADRATIC_CONE n+2 ];
for alpha = alphas
% Objective mu'*x - alpha*s
prob.c = [mu; -alpha];
[rcode,res] = mosekopt('maximize echo(0)',prob,[]);
% Check if the interior point solution is an optimal point
% See https://docs.mosek.com/latest/toolbox/accessing-solution.html about handling solution statuses.
if startsWith(res.rcodestr, 'MSK_RES_ERR') | ~strcmp(res.sol.itr.solsta, 'OPTIMAL')
throw(MException('mosek:solsta', 'Not an optimal solution'));
end
x = res.sol.itr.xx(1:n);
s = res.sol.itr.xx(n+1);
frontier = [frontier; [alpha, mu'*x, sqrt(s)] ];
end
portfolio_3_impact.m
%%
% 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 = 8;
w = 1.0;
mu = [0.07197349 0.15518171 0.17535435 0.0898094 0.42895777 0.39291844 0.32170722 0.18378628]';
x0 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]';
gamma = 0.36;
GT = [ 0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638; ...
0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506; ...
0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914; ...
0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742; ...
0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ; ...
0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187; ...
0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327; ...
0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ];
% Somewhat arbirtrary choice of m
m = 1.0e-2*ones(n,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 echo(0)');
% 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.accs = [ res.symbcon.MSK_DOMAIN_QUADRATIC_CONE 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.accs = [prob.accs repmat([res.symbcon.MSK_DOMAIN_PRIMAL_POWER_CONE, 3, 2, 2.0, 1.0], 1, n) ];
[rcode,res] = mosekopt('maximize echo(0)',prob,[]);
% Check if the interior point solution is an optimal point
% See https://docs.mosek.com/latest/toolbox/accessing-solution.html about handling solution statuses.
if startsWith(res.rcodestr, 'MSK_RES_ERR') | ~strcmp(res.sol.itr.solsta, 'OPTIMAL')
throw(MException('mosek:solsta', 'Not an optimal solution'));
end
x = res.sol.itr.xx(1:n);
t = res.sol.itr.xx(n+(1:n));
portfolio_4_transcost.m
%%
% 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 = 8;
w = 1.0;
mu = [0.07197349 0.15518171 0.17535435 0.0898094 0.42895777 0.39291844 0.32170722 0.18378628]';
x0 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]';
gamma = 0.36;
GT = [ 0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638; ...
0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506; ...
0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914; ...
0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742; ...
0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ; ...
0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187; ...
0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327; ...
0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ];
f = 0.01*ones(n,1);
g = 0.001*ones(n,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 echo(0)');
% 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.accs = [ res.symbcon.MSK_DOMAIN_QUADRATIC_CONE n+1 ];
% Demand y to be integer (hence binary)
prob.ints.sub = 2*n+(1:n);
[rcode,res] = mosekopt('maximize echo(0)',prob,[]);
% Check if the integer solution is an optimal point
% See https://docs.mosek.com/latest/toolbox/accessing-solution.html about handling solution statuses.
if startsWith(res.rcodestr, 'MSK_RES_ERR') | ~strcmp(res.sol.int.solsta, 'INTEGER_OPTIMAL')
throw(MException('mosek:solsta', 'Not an optimal solution'));
end
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
%%
% 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 = 8;
w = 1.0;
mu = [0.07197349 0.15518171 0.17535435 0.0898094 0.42895777 0.39291844 0.32170722 0.18378628]';
x0 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]';
gamma = 0.25;
GT = [ 0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638; ...
0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506; ...
0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914; ...
0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742; ...
0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ; ...
0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187; ...
0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327; ...
0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ];
disp(sprintf('\nMarkowitz portfolio optimization with cardinality constraints'));
for k=1:n
x = MarkowitzWithCardinality(n,mu,GT,x0,w,gamma,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 echo(0)');
% 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.accs = [ res.symbcon.MSK_DOMAIN_QUADRATIC_CONE n+1 ];
% Demand y to be integer (hence binary)
prob.ints.sub = 2*n+(1:n);
[rcode,res] = mosekopt('maximize echo(0)',prob,[]);
% Check if the integer solution is an optimal point
% See https://docs.mosek.com/latest/toolbox/accessing-solution.html about handling solution statuses.
if startsWith(res.rcodestr, 'MSK_RES_ERR') | ~strcmp(res.sol.int.solsta, 'INTEGER_OPTIMAL')
throw(MException('mosek:solsta', 'Not an optimal solution'));
end
x = res.sol.int.xx(1:n);
portfolio_6_factor.m
%%
% File : portfolio_6_factor.m
%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% Description : Implements a basic portfolio optimization model
% of factor type.
%
%%
function portfolio_6_factor()
n = 8;
w = 1.0;
mu = [0.07197349 0.15518171 0.17535435 0.0898094 0.42895777 0.39291844 0.32170722 0.18378628]';
x0 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]';
gammas = [0.24 0.28 0.32 0.36 0.4 0.44 0.48];
B = [ 0.4256, 0.1869; ...
0.2413, 0.3877; ...
0.2235, 0.3697; ...
0.1503, 0.4612; ...
1.5325, -0.2633; ...
1.2741, -0.2613; ...
0.6939, 0.2372; ...
0.5425, 0.2116 ];
S_F = [0.0620, 0.0577; ...
0.0577, 0.0908 ];
theta = [0.0720 0.0508 0.0377 0.0394 0.0663 0.0224 0.0417 0.0459];
%% Compute the small factorization required for the model
P = chol(S_F)';
G_factor = B * P;
G_factor_T = G_factor';
disp(sprintf('\nBasing Markowitz portfolio optimization using a factor model'));
for gamma=gammas
[er, x] = FactorMarkowitz(n,mu,G_factor_T,theta,x0,w,gamma);
disp(sprintf('Expected return: %.4e Std. deviation: %.4e', mu'*x, gamma));
end
%
% Purpose:
% Computes the optimal portfolio for a given risk
%
% Input:
% n: Number of assets
% mu: An n dimmensional vector of expected returns
% G_factor_T: The factor (dense) part of the factorized risk
% theta: specific risk vector
% x0: Initial holdings
% w: Initial cash holding
% gamma: Maximum risk (=std. dev) accepted
%
% Output:
% Optimal expected return and the optimal portfolio
%
function [er, x] = FactorMarkowitz(n,mu,G_factor_T,theta,x0,w,gamma)
[rcode, res] = mosekopt('symbcon echo(0)');
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, G_factor_T*x, sqrt(theta).'x ] in quadratic cone
prob.f = sparse([ zeros(1,n); G_factor_T; diag(sqrt(theta)) ]);
prob.g = [gamma; zeros(size(G_factor_T,1)+n,1)];
prob.accs = [ res.symbcon.MSK_DOMAIN_QUADRATIC_CONE size(G_factor_T,1)+n+1 ];
% Maximize problem and return the objective value
[rcode,res] = mosekopt('maximize echo(0)', prob, []);
% Check if the interior point solution is an optimal point
% See https://docs.mosek.com/latest/toolbox/accessing-solution.html about handling solution statuses.
if startsWith(res.rcodestr, 'MSK_RES_ERR') | ~strcmp(res.sol.itr.solsta, 'OPTIMAL')
throw(MException('mosek:solsta', 'Not an optimal solution'));
end
x = res.sol.itr.xx;
er = mu'*x;
pow1.m
%
% 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.
% Variables number 1,2,3 correspond to x,y,z, variables 4,5 are auxiliary
prob.c = [-1 0 0 1 1];
prob.a = [1 1 0.5 0 0];
prob.blc = [2.0];
prob.buc = [2.0];
prob.blx = [-inf -inf -inf -inf -inf];
prob.bux = [ inf inf inf inf inf];
% Specify the cones as affine conic constraints.
% Two conic constrains with the power cone, both of dimension 3
prob.accs = [res.symbcon.MSK_DOMAIN_PRIMAL_POWER_CONE 3 2 0.2 0.8 res.symbcon.MSK_DOMAIN_PRIMAL_POWER_CONE 3 2 0.4 0.6];
% The matrices such that f * x + g = [x(1), x(2), x(4), x(3), 1, x(5)]
prob.f = sparse( [1, 2, 3, 4, 6], [1, 2, 4, 3, 5], ones(1, 5) );
prob.g = [0 0 0 0 1 0];
% That means
%
% (x(1), x(2), x(4)) \ in PPOW_3(0.2, 0.8)
% (x(3), 1, x(5)) \ in PPOW_3(0.4, 0.6)
%
% which is equivalent to
%
% |x(4)| <= x(1)^0.2 * x(2)^0.8
% |x(5)| <= x(3)^0.4
% Optimize the problem.
[r,res]=mosekopt('maximize',prob);
% Display the primal solution.
res.sol.itr.xx'
qcqo1.m
%
% 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
%
% 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
%
% 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
%
% 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
%%
% 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
%
% 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
%
% 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.accs = [ sym.MSK_DOMAIN_QUADRATIC_CONE 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
%
% 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];
% The scalar variables appear in an affine conic constraint
prob.accs = [res.symbcon.MSK_DOMAIN_QUADRATIC_CONE 3];
prob.f = speye(3);
[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
%
% 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
sdo_lmi.m
%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File : sdo_lmi.m
%
% Purpose : To solve a problem with an LMI and an affine conic constrained problem with a PSD term
%
% minimize Tr [1, 0; 0, 1]*X + x(1) + x(2) + 1
%
% subject to Tr [0, 1; 1, 0]*X - x(1) - x(2) >= 0
% x(1) [0, 1; 1, 3] + x(2) [3, 1; 1, 0] - [1, 0; 0, 1] >> 0
% X >> 0
%%
function sdo_lmi()
[r, res] = mosekopt('symbcon');
symbcon = res.symbcon;
% The scalar part, as in linear optimization examples
prob.c = [1.0 1.0];
prob.cfix = 1.0;
prob.a = sparse([], [], [], 0, 2); % 0 constraints, 2 scalar variables
prob.blc = []; % Bounds
prob.buc = [];
prob.blx = [-inf, -inf];
prob.bux = [inf, inf];
prob.f = sparse([1, 1, 2, 3, 3, 4], ...
[1, 2, 2, 1, 2, 1], ...
[-1, -1, 3, sqrt(2), sqrt(2), 3], ...
4, 2);
prob.g = [0, -1, 0, -1];
% Specify the affine conic structure
prob.accs = [symbcon.MSK_DOMAIN_RPLUS 1 symbcon.MSK_DOMAIN_SVEC_PSD_CONE 3];
% Dimensions of PSD variables
prob.bardim = [2];
% Block triplet format specifying the lower triangular part
% of the symmetric coefficient matrix 'barc':
prob.barc.subj = [1, 1, 1];
prob.barc.subk = [1, 2, 2];
prob.barc.subl = [1, 1, 2];
prob.barc.val = [1, 0, 1];
% Block triplet format specifying the lower triangular part
% of the symmetric coefficient matrix 'barF' for the ACC:
prob.barf.subi = [1, 1];
prob.barf.subj = [1, 1];
prob.barf.subk = [1, 2];
prob.barf.subl = [1, 1];
prob.barf.val = [0, 1];
% Solve with log output
[r, res] = mosekopt('minimize', prob);
% Print the solution
X = zeros(2);
X([1,2,4]) = res.sol.itr.barx;
X = X + tril(X,-1)';
x = res.sol.itr.xx;
X
x
sensitivity.m
%
% 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
%
% 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
%
% 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
%
% 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
