15 List of examples¶
List of examples shipped in the distribution of API for MATLAB:
File |
Description |
|---|---|
A simple problem with one affine conic constraint (ACC) |
|
A simple problem using affine conic constraints |
|
A simple problem using affine conic constraints |
|
A simple conic exponential problem |
|
A simple conic quadratic problem |
|
Solving Stigler’s Nutrition model |
|
A simple problem with disjunctive constraints (DJC) |
|
Demonstrates a small one-facility location problem (CQO) |
|
A simple geometric program (GP) in conic form |
|
A Hello World example |
|
A simple linear problem using msklpopt |
|
A simple linear problem solved with the simplex solver |
|
A simple linear problem using mosekopt |
|
A simple linear problem solved with the simplex solver |
|
A simple mixed-integer conic problem |
|
A simple mixed-integer linear problem |
|
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 |
|
Demonstrate how to modify and re-optimize a linear problem |
|
Demonstrates proper response handling |
Additional examples can be found on the MOSEK website and in other MOSEK publications.
acc1.m
%%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File : acc1.m
%
% Purpose : Tutorial example for affine conic constraints.
% Models the problem:
%
% maximize c^T x
% subject to sum(x) = 1
% gamma >= |Gx+h|_2
%%
function [xx,prosta,solsta] = acc1()
% Example data
c = [2 3 -1]';
gamma = 0.03;
G = [1.5 0.1 0; ...
0.3 0 2.1];
h = [0 0.1]';
n = size(c, 1);
k = size(h, 1);
% Initialize the model
model = mosekmodel(name = "acc1",...
numvar = n);
% Set objective vector
model.objective("maximize", c);
% The constraint sum(x) = 1
model.appendcons(F = ones(1,n), domain = mosekdomain("equal", rhs=1));
% The conic quadratic constraint Fx+g \in \Quad with k+1 rows
model.appendcons(F = sparse([zeros(1,n); G]), ...
g = [gamma; h], ...
domain = mosekdomain("quadratic cone", dim=k+1));
% Solve the problem
model.solve();
% Check if solution is available
[hassol, prosta, solsta] = model.hassolution("interior");
if hassol && solsta == "OPTIMAL"
% Get primal solution
xx = model.getsolution("interior", "x");
% Get dual solution
y = model.getsolution("interior", "y");
disp("Primal solution");
disp(xx);
end
end
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()
% Variables [x1; x2; t1; t2]
model = mosekmodel(name = "affco1", numvar = 4, ...
objsense = "maximize", objective = [0, 0, 1, 1], ...
F = sparse([zeros(1,4); speye(2) zeros(2,2)]), ... % The quadratic cone constraint
g = [1 -0.5 -0.6]', ...
domain = mosekdomain("qcone", dim = 3));
% The power cones added as one block:
model.appendcons(name="pow", ...
F = sparse([1,3,4,4,6], [1,3,1,2,4], ones(1,5)), ...
g = [0 1 0 0.1 1 0]', ...
domain = [mosekdomain("pow", dim = 3, alpha = [1 2]'), ... % Exponents [ 1/3, 2/3 ]
mosekdomain("pow", dim = 3, alpha = 0.25) ]); % Exponents [ 0.25, 0.75 ]
% Linear inequality x_1 - x_2 <= 1
model.appendcons(F = [1 -1 0 0], domain = mosekdomain("less than", rhs = 1));
% Solve and get solution
model.solve();
if model.hassolution("interior")
[xx, prosta, solsta] = model.getsolution("interior", "x");
x = xx(1:2);
t = xx(3:4);
disp(x);
end
end
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 = 3;
z = [3.3, 2.2, 1.3]';
a = [-0.08, -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)
% Variables [t, y1, ..., yn]
model = mosekmodel(...
name = "affco2", ...
numvar = n + 1, ...
objective = [1 zeros(1,n)], ...
objsense = "minimize");
model.varname([1], ["t"]);
% Quadratic cone
model.appendcons(F = diag([0; z]), ...
g = [d; zeros(n,1)], ...
domain = mosekdomain("qcone", dim = n + 1));
% All exponential cones (their number is n)
FExp = sparse([1:3:3*n 3:3:3*n], ...
[2:n+1 ones(1,n)], ...
[ones(1,n) a']);
gExp = repmat([0; 1; 0], n, 1);
model.appendcons(F = FExp, g = gExp, ...
domain = mosekdomain("exp", dim = 3, n = n));
% Solve and get solution
model.solve();
if model.hassolution("interior")
[x, prosta, solsta] = model.getsolution("interior", "x");
t = x(1);
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 API.
%
function [xx,prosta,solsta] = ceo1()
model = mosekmodel(...
name = "ceo1",...
numvar = 3);
% Minimize x1+x2
model.objective("minimize", [1 1 0]);
% Subject to: sum(x) == 1
model.appendcons(F = [1 1 1], domain = mosekdomain("equal", rhs = 1.0));
% And subject to: (x1, x2, x3) \in ExpCone
model.appendcons(F = speye(3), ...
domain = mosekdomain("exp", dim = 3));
% Solve the model and retrieve the solution
model.solve();
if model.hassolution("interior")
[xx,prosta,solsta] = model.getsolution("interior","x");
disp(xx);
end
end
cqo1.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: cqo1.m
%
% Purpose: Demonstrates how to solve the problem
%
% minimize y1 + y2 + y3
% such that
% x1 + x2 + 2.0 x3 = 1.0
% x1,x2,x3 >= 0.0
% and
% (y1,x1,x2) in C_3,
% (y2,y3,x3) in K_3
%
% where C_3 and K_3 are respectively the quadratic and
% rotated quadratic cone of size 3 defined as
% C_3 = { z1,z2,z3 : z1 >= sqrt(z2^2 + z3^2) }
% K_3 = { z1,z2,z3 : 2 z1 z2 >= z3^2 }
%%
function [xx,prosta,solsta] = cqo1()
model = mosekmodel(...
name = "cqo1", ...
numvar = 6);
% Variable is [x, y]
model.objective("minimize", [0 0 0 1 1 1]');
% Linear constraint
model.appendcons(F = [1 1 2 0 0 0], domain = mosekdomain("eq", rhs=1.0));
% Bounds on x
model.appendcons(F = sparse([eye(3) zeros(3)]), domain = mosekdomain("rplus", dim=3));
% Quadratic cone
model.appendcons(F = sparse([1,2,3],[4,1,2],[1,1,1]), domain = mosekdomain("qcone",dim=3));
% Rotated quadratic cone
model.appendcons(F = sparse([1,2,3],[5,6,3],[1,1,1]), domain = mosekdomain("rqcone",dim=3));
model.solve();
if model.hassolution("interior")
[xx,prosta,solsta] = model.getsolution("interior","x");
fprintf("Solution x1,x2,x3: [ %s ]\n", sprintf("%.2g ", xx(1:3)));
fprintf("Solution y1,y2,y3: [ %s ]\n", sprintf("%.2g ", xx(4:6)));
end
end
diet.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: diet.m
%
% Purpose: Solving Stigler's Nutrition model (DIET,SEQ=7)
%
% Source: GAMS Model library,
% Dantzig, G B, Chapter 27.1. In Linear Programming and Extensions.
% Princeton University Press, Princeton, New Jersey, 1963.
%
% Given a set of nutrients, foods with associated nutrient values, allowance of
% nutrients per day, the model will find the cheapest combination of foods
% which will provide the necessary nutrients.
%%
function [xx,prosta,solsta] = diet()
nutrient_unit = [ ...
"thousands", "grams", "grams", ...
"milligrams", "thousand ius", "milligrams", ...
"milligrams", "milligrams", "milligrams" ];
nutrients = [ ...
"calorie", "protein", "calcium", ...
"iron", "vitamin a", "vitamin b1", ...
"vitamin b2", "niacin", "vitamin c" ];
foods = [ ...
"wheat", "cornmeal", "cannedmilk", "margarine", "cheese", ...
"peanut butter", "lard", "liver", "porkroast", "salmon", ...
"greenbeans", "cabbage", "onions", "potatoes", "spinach", ...
"sweet potatos", "peaches", "prunes", "limabeans", "navybeans" ];
daily_allowance = [ ...
3. 70. .8 12. 5. 1.8 2.7 18. 75. ];
nutritive_value = [...
% calorie calcium vitamin a vitamin b2 vitamin c
% protein iron vitamin b1 niacin
44.7 1411 2.0 365 0 55.4 33.3 441 0; % wheat
36 897 1.7 99 30.9 17.4 7.9 106 0; % cornmeal
8.4 422 15.1 9 26 3 23.5 11 60; % cannedmilk
20.6 17 .6 6 55.8 .2 0 0 0; % margarine
7.4 448 16.4 19 28.1 .8 10.3 4 0; % cheese
15.7 661 1 48 0 9.6 8.1 471 0; % peanut butter
41.7 0 0 0 0.2 0 .5 5 0; % lard
2.2 333 .2 139 169.2 6.4 50.8 316 525; % liver
4.4 249 .3 37 0 18.2 3.6 79 0; % porkroast
5.8 705 6.8 45 3.5 1 4.9 209 0; % salmon
2.4 138 3.7 80 69 4.3 5.8 37 862; % greenbeans
2.6 125 4 36 7.2 9 4.5 26 5369; % cabbage
5.8 166 3.8 59 16.6 4.7 5.9 21 1184; % onions
14.3 336 1.8 118 6.7 29.4 7.1 198 2522; % potatoes
1.1 106 0.0 138 918.4 5.7 13.8 33 2755; % spinach
9.6 138 2.7 54 290.7 8.4 5.4 83 1912; % sweet potatos
8.5 87 1.7 173 86.8 1.2 4.3 55 57; % peaches
12.8 99 2.5 154 85.7 3.9 4.3 65 257; % prunes
17.4 1055 3.7 459 5.1 26.9 38.2 93 0; % limabeans
26.9 1691 11.4 792 0 38.4 24.6 217 0] % navybeans
[N,M] = size(nutritive_value);
model = mosekmodel(...
name = "Stinglers Diet Problem", ...
numvar = N);
model.varname([1:N],foods);
model.appendcons(name = "Daily purchase", F = speye(N), domain = mosekdomain("greater than", n=N, rhs=zeros(N,1)));
model.appendcons(name = "Nutrient balance", F = nutritive_value', domain = mosekdomain("greater than", n=M, rhs=daily_allowance'));
model.objective("minimize", ones(N,1));
model.solve();
[ok,prosta,solsta] = model.hassolution("interior");
if ok
fprintf("Solution status: %s, Problem status: %s\n",solsta,prosta);
[xx,prosta,solsta] = model.getsolution("interior","x");
end
end
djc1.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: djc1.m
%
% Purpose: Demonstrates how to solve the problem with two disjunctions:
%
% minimize 2x0 + x1 + 3x2 + x3
% subject to x0 + x1 + x2 + x3 >= -10
% (x0-2x1<=-1 and x2=x3=0) or (x2-3x3<=-2 and x1=x2=0)
% x0=2.5 or x1=2.5 or x2=2.5 or x3=2.5
%%
function [xx,prosta,solsta] = djc1()
model = mosekmodel(name = "djc1", ...
objsense = "minimize", ...
objective = [ 2,1,3,1 ]', ...
numvar = 4);
% A disjunction of two clauses
model.appenddisjunction( [ model.clause(F = [1 -2 0 0 ;
0 0 1 0 ;
0 0 0 1], ...
domain=[ mosekdomain("less than", rhs = [-1]), ... % 1st simple term of 1st clause
mosekdomain("equal", dim = 2, rhs = [0 0]') ]), ... % 2nd simple term of 1st clause
model.clause(F = [0 0 1 -3 ;
1 0 0 0 ;
0 1 0 0], ...
domain = [ mosekdomain("less than", rhs = -2), ... % 1st simple term of 2nd clause
mosekdomain("equal", dim = 2, rhs = [0 0]') ]) ], ... % 2nd simple term of 2nd clause
name = "first-djc");
% A disjunction of four clauses
model.appenddisjunction([ model.clause(F = [1 0 0 0], domain = mosekdomain("equal", rhs = 2.5)),...
model.clause(F = [0 1 0 0], domain = mosekdomain("equal", rhs = 2.5)),...
model.clause(F = [0 0 1 0], domain = mosekdomain("equal", rhs = 2.5)),...
model.clause(F = [0 0 0 1], domain = mosekdomain("equal", rhs = 2.5)) ],...
name = "second-djc");
% The standard liner constraint
model.appendcons(name = "C", F = [1 1 1 1], domain = mosekdomain("greater than", rhs = -10));
model.solve();
if model.hassolution("integer")
[xx,prosta,solsta] = model.getsolution("integer","x");
fprintf("Solution : [%s]\n", sprintf("%g ", xx));
end
end
facility_location.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: facility_location.m
%
% Purpose: Demonstrates a small one-facility location problem.
%
% Given 10 customers placed in a grid we wish to place a facility
% somewhere so that the total sum of distances to customers is
% minimized.
%
% The problem is formulated as a conic optimization problem as follows.
% Let f=(fx,fy) be the (unknown) location of the facility, and let
% c_i=(cx_i,cy_i) be the (known) customer locations; then we wish to
% minimize
% sum_i || f - c_i ||
% where
% ||.||
% denotes the euclidian norm.
% This is formulated as
%
% minimize sum(d_i)
% such that d_i ^ 2 > tx_i ^ 2 + ty_i ^ 2, for all i (a.k.a. (d_i,tx_i,ty_i) in C^3_r)
% tx_i = cx_i - fx, for all i
% ty_i = cy_i - fy, for all i
% d_i > 0, for all i
%%
function [xx,prosta,solsta] = facility_location()
customerloc = [12.0 2.0 ; ...
15.0 13.0 ; ...
10.0 8.0 ; ...
0.0 10.0 ; ...
6.0 13.0 ; ...
5.0 8.0 ; ...
10.0 12.0 ; ...
4.0 6.0 ; ...
5.0 2.0 ; ...
1.0 10.0];
N = size(customerloc,1);
model = mosekmodel(...
name="FacilityLocation",...
numVar=2+N);
% variable(1:N) are the distances from the facility to each customerloc
% variable(N+1:N+2) is the location of the facility
% distance per customer
appendcons(model,...
F = sparse([1:N*3], reshape([ [1:N] ; repmat([N+1;N+2],1,N)],N*3,1)', ones(1,N*3)),...
g = reshape([ zeros(1,N) ; -customerloc' ], 3*N,1),...
domain=mosekdomain("qcone",n=N,dim=3));
% Summed distances
objective(model,...
"minimize", ...
[ 0; 0; ones(N,1)] );
solve(model);
if hassolution(model,"interior")
[xx,prosta,solsta] = getsolution(model,"interior","x");
end
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
function [xx,prosta,solsta] = gp1()
Awall = 200
Afloor = 50
alpha = 2
beta = 10
gamma = 2
delta = 10
% A model with variables [x,y,z,u1,u2]
model = mosekmodel(...
name = "gp1",...
numvar = 5);
model.objective("maximize", [ 1 1 1 0 0]);
% u1 + u2 = 1
model.appendcons(F = [ 0 0 0 1 1 ], domain = mosekdomain("equal", rhs=1.0));
% y + z <= log(Afloor)
model.appendcons(F = [ 0 1 1 0 0 ], domain = mosekdomain("less than", rhs=log(Afloor)));
% Two-sided bounds on x-y and z-y
model.appendcons(F = [ 1 -1 0 0 0 ; ...
0 -1 1 0 0 ], domain = mosekdomain("greater than",dim=2, rhs=[log(alpha) log(gamma)]'));
model.appendcons(F = [ 1 -1 0 0 0 ; ...
0 -1 1 0 0 ], domain = mosekdomain("less than", dim=2, rhs=[log(beta) log(delta)]'));
% Conic constraints
model.appendcons(F = sparse([1 3 3],[4 1 2],[1.0 1.0 1.0]), g = [0 1 log(alpha/Awall)]', domain = mosekdomain("exp"));
model.appendcons(F = sparse([1 3 3],[5 1 3],[1.0 1.0 1.0]), g = [0 1 log(alpha/Awall)]', domain = mosekdomain("exp"));
model.solve();
if model.hassolution("interior")
[xx,prosta,solsta] = model.getsolution("interior","x");
x = xx(1:3);
disp(exp(x));
end
end
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.
%
function [xx,prosta,solsta] = helloworld()
model = mosekmodel(...
name = "helloworld", ...
objsense = "min", ...
objective = [ 1, 2 ]', ...
numvar = 2, ...
F = [ 1 0 ; ...
0 1 ; ...
3 -1 ], ...
domain = [ mosekdomain("greater than", rhs = 2), ... % x >= 2, 1st row of F
mosekdomain("greater than", rhs = 3), ... % y >= 3, 2nd row of F
mosekdomain("less than", rhs = 1) ]); % 3x-y <= 1, 3rd row of F
model.solve();
if model.hassolution("interior")
[xx,prosta,solsta] = model.getsolution("any", "x");
xx
end
end
lo1.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: lo1.m
%
% Purpose: Demonstrates how to solve the problem
%
% maximize 3*x0 + 1*x1 + 5*x2 + x3
% such that
% 3*x0 + 1*x1 + 2*x2 = 30,
% 2*x0 + 1*x1 + 3*x2 + 1*x3 > 15,
% 2*x1 + + 3*x3 < 25
% and
% x0,x1,x2,x3 > 0,
% 0 < x1 < 10
%%
function [xx,prosta,solsta] = lo1()
model = mosekmodel(...
name = "lo1", ...
objsense = "maximize", ...
objective = [ 3 1 5 1 ]', ...
numvar = 4, ...
F = [ 3 1 2 0 ; ...
2 1 3 1 ; ...
0 2 0 3 ; ...
1 0 0 0 ; ...
0 1 0 0 ; ...
0 0 1 0 ; ...
0 0 0 1 ; ...
0 1 0 0 ;], ...
domain = [ mosekdomain("equal", rhs=30), ...
mosekdomain("greater than", rhs=15), ...
mosekdomain("less than", rhs=25), ...
mosekdomain("nonnegative", n=4), ...
mosekdomain("less than", rhs=10) ]);
model.varname([1:4],["x","y","z","w"]);
model.solve();
if model.hassolution("interior")
[xx,prosta,solsta] = model.getsolution("interior","x");
end
end
lo1_simplex.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: lo1_simplex.m
%
% Purpose: Demonstrates how to solve the problem
%
% maximize 3*x0 + 1*x1 + 5*x2 + x3
% such that
% 3*x0 + 1*x1 + 2*x2 = 30,
% 2*x0 + 1*x1 + 3*x2 + 1*x3 > 15,
% 2*x1 + + 3*x3 < 25
% and
% x0,x1,x2,x3 > 0,
% 0 < x1 < 10
%
% using the linear/simplex toolbox.
%%
function [xx,prosta,solsta] = lo1_simplex()
model = moseklinmodel(...
name = "lo1", ...
numvar = 6, ...
objsense = "maximize", ...
objective = [ 3 1 5 1 0 0 ]', ...
A = [ 3 1 2 0 0 0 ; ...
2 1 3 1 -1 0 ; ...
0 2 0 3 0 -1 ], ...
b = [ 30 15 25 ]', ...
blx = [ 0.0 0.0 0.0 0.0 0.0 -inf ]', ...
bux = [ inf 10.0 inf inf inf 0.0 ]', ...
varnames = ["x1" "x2" "x3" "x4" "y1" "y2"]',...
conname = ["c1" "c2" "c3"]');
% Choose the simplex optimizer to solve with
model.solve(param = [ "MSK_IPAR_OPTIMIZER", "MSK_OPTIMIZER_FREE_SIMPLEX" ]);
% Check the solution status and fetch the solution
[hassol, prosta, solsta] = model.hassolution();
if hassol && solsta == "OPTIMAL"
fprintf("Solution, problem status = %s, solution status = %s:\n", prosta, solsta);
% Primal variables
result = model.getsolution("x");
x = result(1:4);
s = result(5:6);
fprintf(" x[%d] = %.4f\n",[1:4 ; x']);
fprintf("Slacks:\n");
fprintf(" s[%d] = %.4f\n",[1:2 ; s']);
% Other elements of the solution
fprintf("Duals:\n");
fprintf(" y[%d] = %.4f\n",[1:3 ; model.getsolution("y")']);
fprintf("slx[%d] = %.4f\n",[1:6 ; model.getsolution("slx")']);
fprintf("sux[%d] = %.4f\n",[1:6 ; model.getsolution("sux")']);
fprintf("Status keys of constraints:\n");
fprintf("skc[%s] = %s\n",[1:3 ; model.getsolution("skc")']);
fprintf("Status keys of variables:\n");
fprintf("skx[%s] = %s\n",[1:6 ; model.getsolution("skx")']);
end
end
lo2.m
%%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File : lo2.m
%
% Purpose: Demonstrates how to solve the problem
%
% maximize 3*x0 + 1*x1 + 5*x2 + x3
% such that
% 3*x0 + 1*x1 + 2*x2 = 30,
% 2*x0 + 1*x1 + 3*x2 + 1*x3 > 15,
% 2*x1 + + 3*x3 < 25
% and
% x0,x1,x2,x3 > 0,
% 0 < x1 < 10
%%
function [xx,prosta,solsta] = lo2()
model = mosekmodel(name = "lo2", ...
objsense = "maximize", ...
objective = [ 3 1 5 1 ]', ...
numvar = 4);
model.appendcons(name="con-eq30", F = [ 3 1 2 0 ], domain = mosekdomain("equal", rhs=30));
model.appendcons(name="con-gt14", F = [ 2 1 3 1 ], domain = mosekdomain("greater than", rhs=15));
model.appendcons(name="con-lt25", F = [ 0 2 0 3 ], domain = mosekdomain("less than", rhs=25));
model.appendcons(name="con-nneg", F = speye(4), domain = mosekdomain("nonnegative", n=4));
model.appendcons(name="con-lt10", F = [ 0 1 0 0 ], domain = mosekdomain("less than", rhs="10"));
model.varname([1:4],["x","y","z","w"]);
model.solve();
if model.hassolution("interior")
[xx,prosta,solsta] = model.getsolution("interior","x");
end
end
lo2_simplex.m
%%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File : lo2_simplex.m
%
% Purpose : Demonstrates how to solve small linear
% optimization problem using the MOSEK API for Matlab
% in multiple calls.
%%
function [xx,prosta,solsta] = lo2_simplex()
% Create an empty model
model = moseklinmodel(name = "lo2");
% Add some variables with bounds
model.appendvars(6,...
bl = [ 0.0 0.0 0.0 0.0 0.0 -inf ]', ...
bu = [ inf 10.0 inf inf inf 0.0 ]');
model.varname([1:6],["x1" "x2" "x3" "x4" "y1" "y2"]);
% Set the objective
model.objective("maximize", [ 3 1 5 1 0 0]', objfixterm = 0.0);
% Constraint 3x1 + x2 + 2x3 == 30
model.appendcons([ 3 1 2 0 0 0], [30]);
% The two remaining constraints
model.appendcons([ 2 1 3 1 -1 0 ; ...
0 2 0 3 0 -1 ], ...
[ 15 25 ]');
% Constraint names
model.conname([1:3],["c1" "c2" "c3"]);
% Choose the simplex optimizer
model.solve(param = [ "MSK_IPAR_OPTIMIZER", "MSK_OPTIMIZER_FREE_SIMPLEX" ]);
% Check the solution status and fetch the solution
[hassol, prosta, solsta] = model.hassolution();
if hassol && solsta == "OPTIMAL"
fprintf("Solution, problem status = %s, solution status = %s:\n", prosta, solsta);
% Primal variables
result = model.getsolution("x");
x = result(1:4);
s = result(5:6);
fprintf(" x[%d] = %.4f\n",[1:4 ; x']);
fprintf("Slacks:\n");
fprintf(" s[%d] = %.4f\n",[1:2 ; s']);
% Other elements of the solution
fprintf("Duals:\n");
fprintf(" y[%d] = %.4f\n",[1:3 ; model.getsolution("y")']);
fprintf("slx[%d] = %.4f\n",[1:6 ; model.getsolution("slx")']);
fprintf("sux[%d] = %.4f\n",[1:6 ; model.getsolution("sux")']);
fprintf("Status keys of constraints:\n");
fprintf("skc[%s] = %s\n",[1:3 ; model.getsolution("skc")']);
fprintf("Status keys of variables:\n");
fprintf("skx[%s] = %s\n",[1:6 ; model.getsolution("skx")']);
end
end
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()
%The full variable is [t; x; y]
model = mosekmodel(name = "mi-conic", ...
numvar = 3, objective = [1 0 0], objsense = "min", ...
F = [ speye(3); ...
0 1 0; ...
0 0 0; ...
0 0 1 ], ...
g = [0 0 0 -3.8 1 0], ...
domain = [mosekdomain("qcone", dim = 3), ...
mosekdomain("exp")], ...
intvars = [2, 3]); % Specify the indices of integer variables
% It is as always possible (but not required) to input an initial solution
% to start the mixed-integer solver.
model.setsolution("x", [100, 9, -1]);
model.solve();
x = model.getsolution("integer", "x");
fprintf("Optimal (x,y) = (%.2e, %.2e)\n", x(2), x(3));
end
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 [xx,prosta,solsta] = milo1()
model = mosekmodel(name="milo1", numvar=2, ...
intvars=[1 2]'); % Specify indices of integer variables
model.objective("maximize", [ 1.0 0.64 ]);
model.appendcons(F = speye(2), domain = mosekdomain("rplus",dim=2));
model.appendcons(F = [50.0 31.0], domain = mosekdomain("less than", rhs=250.0));
model.appendcons(F = [4.0 -2.0], domain = mosekdomain("greater than", rhs=-4));
model.solve();
% Access the integer solution
if model.hassolution("integer")
[xx,prosta,solsta] = model.getsolution("integer", "x");
fprintf("Solution: %s\n", sprintf("%g ", xx));
else
disp("Solve failed");
end
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()
model = mosekmodel(name = "init sol", ...
numvar = 4, objective = [7 10 1 5], objsense = "max", ...
F = [1 1 1 1; eye(4)], ...
domain = [mosekdomain("eq", rhs = 2.5), ...
mosekdomain("gt", dim = 4, rhs = zeros(4, 1))], ...
intvars = [1, 2, 3]); % Specify the indices of integer variables
% Specify start guess for the integer variables.
model.setsolution("x", [1 1 0 nan]);
% Request constructing the solution from integer variable values
model.solve(param = ["MSK_IPAR_MIO_CONSTRUCT_SOL", "1"]);
% Display the optimal solution.
if model.hassolution("integer")
x = model.getsolution("integer", "x");
fprintf("x = [%s]\n", sprintf("%g ", x));
% Was the initial solution used ?
fprintf("Construct solution used? %d\n", model.info.MSK_IINF_MIO_CONSTRUCT_SOLUTION)
fprintf("Construct solution objective: %f\n", model.info.MSK_DINF_MIO_CONSTRUCT_SOLUTION_OBJ);
end
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)
n = size(F,2);
k = size(g,1);
m = size(A,1);
% A model with variables x(n) and t(1)
model = mosekmodel(name = "norm_lse", numvar = n + 1, objective = [zeros(n, 1); 1]);
% Linear constraints
model.appendcons(F = [A, zeros(m, 1)], domain = mosekdomain("equal", dim = m, rhs = b));
% Affine conic constraint
model.appendcons(F = [sparse(1,n), 1; F, sparse(k,1)], g = [0; -g], ...
domain = mosekdomain("quadratic cone", dim = k + 1));
model.solve();
xx = model.getsolution("any", "x");
x = 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)
n = size(F,2);
k = size(g,1);
m = size(A,1);
% A model with variables x(n) and t1, t2
model = mosekmodel(name = "norm_lse_reg", numvar = n + 2, objective = [zeros(n, 1); 1; gamma]);
% Linear constraints
model.appendcons(F = [A, zeros(m, 2)], domain = mosekdomain("equal", dim = m, rhs = b));
% Affine conic constraint for \|Fx-g\|_2
model.appendcons(F = [sparse(1,n), 1, 0; ...
F, sparse(k,2)], ...
g = [0; -g], ...
domain = mosekdomain("quadratic cone", dim = k + 1));
% Affine conic constraint for \|x\|_2
model.appendcons(F = [sparse(1, n+1), 1; ...
speye(n), sparse(n,2) ], ...
domain = mosekdomain("quadratic cone", dim = n + 1));
model.solve();
xx = model.getsolution("any", "x");
x = 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)
n = size(F,2);
k = size(g,1);
m = size(A,1);
% A model with variables x(n) and t1, t2
model = mosekmodel(name = "norm_lse_reg_quad", numvar = n + 2, objective = [zeros(n, 1); 1; gamma]);
% Linear constraints
model.appendcons(F = [A, zeros(m, 2)], domain = mosekdomain("equal", dim = m, rhs = b));
% Affine conic constraint for \|Fx-g\|_2^2
model.appendcons(F = [sparse(1,n) 1 0; ...
sparse(1, n+2); ...
F sparse(k, 2)], ...
g = [0; 0.5; -g], ...
domain = mosekdomain("rotated quadratic cone", dim = k + 2));
% Affine conic constraint for \|x\|_2^2
model.appendcons(F = [sparse(1, n+1) 1; ...
sparse(1, n+2); ...
speye(n) sparse(n, 2) ], ...
g = [0; 0.5; zeros(n, 1)], ...
domain = mosekdomain("rotated quadratic cone", dim = n + 2));
model.solve();
xx = model.getsolution("any", "x");
x = 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)
n = size(F,2);
k = size(g,1);
m = size(A,1);
% A model with variables x(n), u(n) and t1, t2
model = mosekmodel(name = "norm_lse_lasso", numvar = n + n + 2, objective = [zeros(2*n, 1); 1; gamma]);
% Linear constraints
model.appendcons(F = [A, zeros(m, 2)], domain = mosekdomain("equal", dim = m, rhs = b));
% u >=. |x|
model.appendcons(F = [speye(n) speye(n) sparse(n, 2); ...
-speye(n) speye(n) sparse(n, 2)], ...
domain = mosekdomain("greater than", dim = 2*n, rhs = zeros(2*n, 1)));
% t2 >= sum(u)
model.appendcons(F = [sparse(1, n), -ones(1, n), 0, 1], ...
domain = mosekdomain("greater than", rhs = 0));
% Affine conic constraint for \|Fx-g\|_2
model.appendcons(F = [sparse(1, 2*n) 1 0; ...
F sparse(k, n + 2)], ...
g = [0; -g], ...
domain = mosekdomain("quadratic cone", dim = k + 1));
model.solve();
xx = model.getsolution("any", "x");
x = xx(1:n);
end
% P-norm minimization
% minimize \|Fx-g\|_p
function x = norm_p_norm(F,g,A,b,p)
n = size(F,2);
k = size(g,1);
m = size(A,1);
% A model with variables x(n), r(k) and t(1)
model = mosekmodel(name = "norm_p_norm", numvar = n + k + 1, objective = [zeros(n + k, 1); 1]);
% Linear constraints
model.appendcons(F = [A, sparse(m, k + 1)], domain = mosekdomain("equal", dim = m, rhs = b));
% t == sum(r)
model.appendcons(F = [sparse(1, n), ones(1, k), -1], ...
domain = mosekdomain("equal", rhs = 0));
% 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
Fcon = M * [sparse(k, n), speye(k), sparse(k,1); ...
sparse(k, n+k), ones(k, 1); ...
F, sparse(k, k+1)];
gcon = M * [sparse(2*k, 1); -g];
model.appendcons(F = Fcon, g = gcon, ...
domain = mosekdomain("power cone", n = k, dim = 3, alpha = [1; p-1]));
model.solve();
xx = model.getsolution("any", "x");
x = 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(url, cert)
% Here we can set up a model
model = mosekmodel(...
objsense = "min", ...
objective = [ 1, 2 ]', ...
numvar = 2, ...
F = [ 1 0 ; ...
0 1 ; ...
3 -1 ], ...
domain = [ mosekdomain("greater than", rhs = 2), ...
mosekdomain("greater than", rhs = 3), ...
mosekdomain("less than", rhs = 1) ]);
% Set up the certificate path, if using TLS, otherwise ignore
if exist('cert','var')
param = ["MSK_SPAR_REMOTE_TLS_CERT_PATH", cert ];
else
param = [];
end
% Optimize using the remote server
model.solve(param = param, optserver = url);
% Use the optimal solution
if model.hassolution("interior")
xx = model.getsolution("interior", "x");
disp(xx);
end
parameters.m
%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File : parameters.m
%
% Purpose: Demonstrates setting solver parameters and
% retrieving information items.
%
function parameters()
% Set up a small problem as an example
model = mosekmodel(objsense = "min", ...
objective = [ 1, 2 ]', ...
numvar = 2, ...
F = [ 1 0 ; ...
0 1 ; ...
3 -1 ], ...
domain = [ mosekdomain("greater than", rhs = 2), ...
mosekdomain("greater than", rhs = 3), ...
mosekdomain("less than", rhs = 1) ]);
% Solve with a list of parameters
model.solve(param = ["MSK_IPAR_LOG", "1", ... % Set log level (integer parameter)
"MSK_IPAR_CACHE_LICENSE", "MSK_OFF", ... % Do not keep the license (integer parameter)
"MSK_DPAR_INTPNT_CO_TOL_REL_GAP", "1.0e-7" ]); % Set relative gap tolerance (double parameter)
% Demonstrate information items after optimization
fprintf("Optimizer time %.3f\n", model.info.MSK_DINF_OPTIMIZER_TIME);
fprintf("#iterations %d\n", model.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
[model] = testProblem();
% Perform the optimization.
% For comparison with the following code, we will also
% ask MOSEK to print the infeasibility certificate to the log
% following optimization.
model.solve(param = ["MSK_IPAR_INFEAS_REPORT_AUTO", "MSK_ON"]);
% Check problem status
[hassol, prosta, solsta] = model.hassolution("interior");
if hassol && prosta == "PRIM_INFEAS"
% Set the tolerance at which we consider a dual value as essential
eps = 1e-7;
% Obtain the dual values (contianing certificate)
y = model.getsolution("interior", "y");
% List all certificate entries with (sufficiently) nonzero dual values
disp("Constraint rows participating in infeasibility: ");
cert = find(abs(y) > eps);
disp(cert);
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 [model] = testProblem()
model = mosekmodel(name = "pinfeas", numvar = 7, ...
objsense = "minimize", objective = [1 2 5 2 1 2 1]');
model.appendcons(name = "s", ...
F = sparse([1,1,2,2,3,3,3], [1,2,3,4,5,6,7], [1,1,1,1,1,1,1]), ...
domain = mosekdomain("less than", dim = 3, rhs=[200, 1000, 1000]));
model.appendcons(name = "d", ...
F = sparse([1,1,2,3,3,4,4], [1,5,2,3,6,4,7], [1,1,1,1,1,1,1]), ...
domain = mosekdomain("equals", dim = 4, rhs=[1100, 200, 500, 500]));
model.appendcons(name = "low", F = speye(7), domain = mosekdomain("rplus", dim = 7));
portfolio_1_basic.m
%
% File : portfolio_1_basic.m
%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% Purpose : Implements a basic portfolio optimization model.
%
function [xx,prosta,solsta] = portfolio_1_basic()
% 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
n = 8;
w = 59;
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;
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 ];
m = size(GT,1);
model = mosekmodel(...
name="Basic Markowitz", ...
numvar = n, ...
objsense = "maximize", ...
objective = mu);
% Bound on the variables (holdings). Shortselling is not allowed.
model.appendcons(F = eye(n), domain = mosekdomain("nonnegative", n = n));
% The amount invested must be identical to initial wealth
model.appendcons(F = ones(1, n), domain = mosekdomain("equal", rhs=w+sum(x0)));
% Imposes a bound on the risk
model.appendcons(F = [zeros(1,n); GT], g = [gamma; zeros(m,1)],...
domain = mosekdomain("qcone", dim = m+1));
model.solve();
% Check if the interior point solution is an optimal point
% See https://docs.mosek.com/latest/matlabapi/accessing-solution.html about handling solution statuses.
[ok,prosta,solsta] = model.hassolution("interior");
if ~ ok || solsta ~= "OPTIMAL"
error("Unexpected solution status");
else
xx = model.getsolution("interior","x");
disp("\n-----------------------------------------------------------------------------------");
disp("Basic Markowitz portfolio optimization");
disp("-----------------------------------------------------------------------------------\n");
fprintf("Expected return: %.4e Std. deviation: %.4e\n", mu' * xx, gamma);
end
end
portfolio_2_frontier.m
%
% File : portfolio_2_frontier.m
%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% Purpose : Implements a basic portfolio optimization model.
% Computes points on the efficient frontier.
%
function [xx,prosta,solsta] = portfolio_2_frontier()
% 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)
n = 8;
w = 1.0;
mu = [0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379]';
x0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.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 ]';
niter = size(alphas,1);
frontier_mux = zeros(niter,1);
frontier_s = zeros(niter,1);
for i = 1:niter
alpha = alphas(i);
model = mosekmodel(name = "Efficient frontier", numvar = n+1);
% Defines the variables (holdings). Shortselling is not allowed.
model.appendcons(F = speye(n), domain = mosekdomain("greater than", dim = n, rhs = zeros(n,1)));
x_idxs = [1:n];
s_idxs = n+1;
model.appendcons(name = "budget", F = [ones(1,n) 0], domain = mosekdomain("equal", rhs = w + sum(x0)));
% Computes the risk
model.appendcons(name="variance", ...
F=[ sparse(1,n+1,[1.0],[1],[s_idxs]) ; ...
sparse(1,n+1) ; ...
GT zeros(n,1)], ...
g=[0; 0.5; zeros(n,1)], ...
domain = mosekdomain("rqcone", dim = n+2));
% Define objective as a weighted combination of return and variance
model.objective("maximize", [ mu ; -alpha ]);
% Solve the problem for the current alpha
model.solve();
% Check if the solution is an optimal point
% See https://docs.mosek.com/latest/matlabapi/accessing-solution.html about handling solution statuses.
[ok,prosta,solsta] = model.hassolution("interior");
if ~ ok || solsta ~= "OPTIMAL"
error("Unexpected solution status");
else
xx = model.getsolution("interior","x");
x = xx(x_idxs,1);
s = xx(s_idxs,1);
frontier_mux(i,1) = mu' * x;
frontier_s(i,1) = s;
end
end
fprintf("\n-----------------------------------------------------------------------------------\n");
fprintf("Efficient frontier") ;
fprintf("\n-----------------------------------------------------------------------------------\n");
fprintf("\t%-12s %-12s %-12s\n", "alpha", "return", "std. dev.");
sqrt_frontier_s = sqrt(frontier_s);
for i = 1:niter
fprintf("\t%-12.4f %-12.4e %-12.4e\n", alphas(i,1), frontier_mux(i,1), frontier_s(i,1));
end
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 [x,prosta,solsta] = portfolio_3_impact()
%
% 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 portfolio (x)
%
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);
% A model with variables (x, t)
model = mosekmodel(name = "Impact", ...
numvar = n+n, ...
objsense = "maximize", ...
objective = [mu; zeros(n,1)]);
% Linear constraint
% [ e' m' ] * [ x; t ] = w + e'*x0
model.appendcons(name = "budget", F = [ ones(1,n), m' ], ...
domain = mosekdomain("equal", rhs = w + sum(x0)));
% No shortselling, x >= 0 and no other bounds
model.appendcons(F = speye(n, 2*n), ...
domain = mosekdomain("greater than", dim = n, rhs = zeros(n,1)));
% Affine conic constraint representing [ gamma, GT*x ] in quadratic cone
model.appendcons(F = [sparse(1, 2*n); ...
GT, sparse(n, n)], ...
g = [gamma; zeros(n, 1)], ...
domain = mosekdomain("qcone", dim = n + 1));
% Power cone constraints 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
model.appendcons(name = "impact", ...
F = sparse(fi, fj, fv, 3*n, 2*n), g = g, ...
domain = mosekdomain("power cone", dim = 3, n = n, alpha = [2.0, 1.0]));
model.solve();
% Check if the solution is an optimal point
% See https://docs.mosek.com/latest/matlabapi/accessing-solution.html about handling solution statuses.
[ok,prosta,solsta] = model.hassolution("interior");
if ~ ok || solsta ~= "OPTIMAL"
error("Unexpected solution status");
else
xx = model.getsolution("interior","x");
x = xx(1:n);
end
end
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 [x,prosta,solsta] = portfolio_4_transcost()
%
% 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 portfolio
%
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);
% Upper bound on the traded amount
u = w + sum(x0);
% A model with variables (x, z, y)
model = mosekmodel(name = "Impact", ...
numvar = n+n+n, ...
objsense = "maximize", ...
objective = [mu; zeros(n,1); zeros(n,1)], ...
intvars = 2*n+1:3*n); % The variables y are integer
% y is binary, ie. in [0,1]
model.appendcons(F = [sparse(n,2*n) speye(n,n)], domain = mosekdomain("gt", dim = n, rhs = zeros(n,1)));
model.appendcons(F = [sparse(n,2*n) speye(n,n)], domain = mosekdomain("lt", dim = n, rhs = ones(n,1)));
% Linear constraints
% [ e' g' f' ] [ x ] = w + e'*x0
% [ I -I 0 ] * [ z ] <= x0
% [ I I 0 ] [ y ] >= x0
% [ 0 I -U ] <= 0
% Linear budget constraint
model.appendcons(name = "budget", F = [ ones(1,n), g', f' ], ...
domain = mosekdomain("equal", rhs = u));
% Linear absolute value model z >= |x-x0|
model.appendcons(name = "abs1", F = [ speye(n) -speye(n) sparse(n,n) ], ...
domain = mosekdomain("lt", dim = n, rhs = x0));
model.appendcons(name = "abs2", F = [ speye(n) speye(n) sparse(n,n) ], ...
domain = mosekdomain("gt", dim = n, rhs = x0));
% "big-M" constraint z<=Uy
model.appendcons(name = "bigM", F = [ sparse(n,n) speye(n) -u*speye(n) ], ...
domain = mosekdomain("lt", dim = n, rhs = zeros(n,1)));
% No shortselling, x >= 0
model.appendcons(F = speye(n, 3*n), ...
domain = mosekdomain("greater than", dim = n, rhs = zeros(n,1)));
% Affine conic constraint representing [ gamma, GT*x ] in quadratic cone
model.appendcons(F = [sparse(1, 3*n); ...
GT, sparse(n, 2*n)], ...
g = [gamma; zeros(n, 1)], ...
domain = mosekdomain("qcone", dim = n + 1));
model.solve();
% Check if the solution is an optimal point
% See https://docs.mosek.com/latest/matlabapi/accessing-solution.html about handling solution statuses.
[ok,prosta,solsta] = model.hassolution("integer");
if ~ ok || solsta ~= "INTEGER_OPTIMAL"
error("Unexpected solution status");
else
xx = model.getsolution("integer","x");
x = xx(1:n);
y = xx(2*n+(1:n));
z = xx(n+(1:n));
end
fprintf("\nMarkowitz portfolio optimization with transactions cost")
fprintf("Expected return: %.4e Std. deviation: %.4e Transactions cost: %.4e", ...
mu'*x, gamma, f'*y+g'*z)
end
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 ];
%
% 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 portfolio
%
fprintf("\nMarkowitz portfolio optimization with cardinality constraints\n");
for k=1:n
% Upper bound on the traded amount
u = w + sum(x0);
% A model with variables (x, z, y)
model = mosekmodel(name = "Transaction costs", ...
numvar = n+n+n, ...
objsense = "maximize", ...
objective = [mu; zeros(n,1); zeros(n,1)], ...
intvars = 2*n+1:3*n); % The variables y are integer
% 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
% Linear budget constraint
model.appendcons(name = "budget", F = [ ones(1,n), zeros(1, 2*n) ], ...
domain = mosekdomain("equal", rhs = u));
% Linear absolute value model z >= |x-x0|
model.appendcons(name = "abs1", F = [ speye(n) -speye(n) sparse(n,n) ], ...
domain = mosekdomain("lt", dim = n, rhs = x0));
model.appendcons(name = "abs2", F = [ speye(n) speye(n) sparse(n,n) ], ...
domain = mosekdomain("gt", dim = n, rhs = x0));
% "big-M" constraint z<=Uy
model.appendcons(name = "bigM", F = [ sparse(n,n) speye(n) -u*speye(n) ], ...
domain = mosekdomain("lt", dim = n, rhs = zeros(n,1)));
% Cardinality bound sum(y) <= k
model.appendcons(name = "card", F = [ zeros(1, 2*n) ones(1, n) ], ...
domain = mosekdomain("lt", rhs = k));
% No shortselling, x >= 0
model.appendcons(F = speye(n, 3*n), ...
domain = mosekdomain("greater than", dim = n, rhs = zeros(n,1)));
% Affine conic constraint representing [ gamma, GT*x ] in quadratic cone
model.appendcons(F = [sparse(1, 3*n); ...
GT, sparse(n, 2*n)], ...
g = [gamma; zeros(n, 1)], ...
domain = mosekdomain("qcone", dim = n + 1));
model.solve(quiet = true);
% Check if the solution is an optimal point
% See https://docs.mosek.com/latest/matlabapi/accessing-solution.html about handling solution statuses.
[ok,prosta,solsta] = model.hassolution("integer");
if ~ ok || solsta ~= "INTEGER_OPTIMAL"
error("Unexpected solution status");
else
xx = model.getsolution("integer","x");
x = xx(1:n);
fprintf("Bound: %d Expected return: %.4e Solution: [ %s ]\n", k, mu'*x, sprintf("%.2g ", x));
end
end
end
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()
%
% 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
%
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';
fprintf("Basing Markowitz portfolio optimization using a factor model");
for gamma = gammas
% A model with variable x and maximized return
model = mosekmodel(name = "Factor risk", ...
numvar = n, ...
objsense = "maximize", ...
objective = mu);
% Linear budget constraint
model.appendcons(name = "budget", F = ones(1,n), ...
domain = mosekdomain("equal", rhs = w + sum(x0)));
% No shortselling, x >= 0
model.appendcons(F = speye(n), ...
domain = mosekdomain("greater than", dim = n, rhs = zeros(n,1)));
% An affine conic constraint: [gamma, G_factor_T*x, sqrt(theta).'x ] in quadratic cone
model.appendcons(name = "risk", ...
F = sparse([ zeros(1,n); G_factor_T; diag(sqrt(theta)) ]), ...
g = [gamma; zeros(size(G_factor_T, 1) + n, 1)], ...
domain = mosekdomain("qcone", dim = size(G_factor_T, 1) + n + 1 ));
model.solve(quiet = true);
% Check if the solution is an optimal point
% See https://docs.mosek.com/latest/matlabapi/accessing-solution.html about handling solution statuses.
[ok,prosta,solsta] = model.hassolution("interior");
if ~ ok || solsta ~= "OPTIMAL"
error("Unexpected solution status");
else
x = model.getsolution("interior","x");
er = mu'*x;
fprintf("Expected return: %.4e Std. deviation: %.4e\n", er, gamma);
end
end
end
pow1.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: pow1.m
%
% Purpose: Demonstrates 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 [xx,prosta,solsta] = pow1()
model = mosekmodel(...
name = "pow1",...
numvar = 5);
model.objective("maximize", [-1 0 0 1 1]);
model.appendcons(F = [1 1 0.5], domain = mosekdomain("equal",rhs=2));
model.appendcons(F = sparse([1 2 3],[1 2 4],[1.0 1.0 1.0]), domain = mosekdomain("pow", dim=3, alpha=[0.2 0.8]));
model.appendcons(F = sparse([1 3],[3 5],[1.0 1.0]), g = [0 1 0]', domain = mosekdomain("pow", dim=3, alpha=[0.4 0.6]));
model.solve();
if model.hassolution("interior")
[xx,prosta,solsta] = model.getsolution("interior","x");
end
end
reoptimization.m
%%
% Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File : reoptimization.m
%
% Purpose: Demonstrates how to solve a linear
% optimization problem using the Conic Toolbox for Matlab
% and modify and re-optimize the problem.
%%
function reoptimization()
model = moseklinmodel(...
name = "reoptimization",...
objsense = "maximize", ...
numvar = 7, ...
blx = [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ]', ...
bux = [ inf inf inf inf inf inf inf ]', ...
b = [ 100000.0 50000.0 60000.0 ]', ...
... % x1 x2 x3 s1 s2 s3 s4
objective = [ 1.5 2.5 3.0 0.0 0.0 0.0 0.0 ]', ...
A = [ 2.0 4.0 3.0 1.0 0.0 0.0 0.0 ; ...
3.0 2.0 3.0 0.0 1.0 0.0 0.0 ; ...
2.0 3.0 2.0 0.0 0.0 1.0 0.0 ]);
model.solve();
[result,prosta,solsta] = model.getsolution("x"); x = result';
display(x);
% Alter A and reoptimize
fprintf("Update first row non-zeros\n");
model.setrows([ 3.0 4.0 3.0 0.0 0.0 0.0 0.0 ], first=1);
model.solve();
[result,prosta,solsta] = model.getsolution("x"); x = result';
display(x);
% Add variable
fprintf("Add variable\n");
model.appendvars(1, ...
bl = [ 0.0 ]', ...
A = [ 4.0 0.0 1.0 ]', ...
c = [ 1.0 ]);
model.solve(param = [ "MSK_IPAR_OPTIMIZER","MSK_OPTIMIZER_FREE_SIMPLEX" ]);
[result,prosta,solsta] = model.getsolution("x"); x = result';
display(x);
% Add constraint
fprintf("Add constraints\n");
% x1 x2 x3 s1 s2 s3 s4 x4
model.appendcons([ 1.0 2.0 1.0 0.0 0.0 0.0 1.0 1.0 ], ...
[ 30000.0 ]);
model.solve(param = [ "MSK_IPAR_OPTIMIZER","MSK_OPTIMIZER_FREE_SIMPLEX" ]);
[result,prosta,solsta] = model.getsolution("x"); x = result';
display(x);
% Change bounds
fprintf("Change bounds\n");
model.setb([ 80000.0 40000.0 50000.0 22000.0 ]');
model.solve(param = [ "MSK_IPAR_OPTIMIZER","MSK_OPTIMIZER_FREE_SIMPLEX" ]);
[result,prosta,solsta] = model.getsolution("x"); x = result';
display(x);
end
response.m
%%
% Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
%
% File: response.m
%
% Purpose: This example demonstrates the workflow for proper response handling
%
function response(args)
% Ceate a simple model for the sake of the example
model = mosekmodel(objsense = "min", ...
objective = [ 1, 2 ]', ...
numvar = 2, ...
F = [ 1 0 ; ...
0 1 ; ...
3 -1 ], ...
domain = [ mosekdomain("greater than", rhs = 2), ...
mosekdomain("greater than", rhs = 3), ...
mosekdomain("less than", rhs = 1) ]);
% Solve the model
try
model.solve();
catch ME
warning("An error during optimization; handle it here.");
rethrow(ME);
end
% We check if the interior-point solution exists and what status it has
[exists, prosta, solsta] = model.hassolution("interior");
if exists
disp("Solved the problem with statuses:");
disp(prosta);
disp(solsta);
switch solsta
case "OPTIMAL"
disp("Optimal solution found:");
x = model.getsolution("interior");
disp(x);
case "PRIM_INFEAS_CER"
disp("The problem is primal infeasible.");
case "DUAL_INFEAS_CER"
disp("The problem is dual infeasible.");
case "UNKNOWN"
disp("Solution status UNKNOWN. This could indicate numerical issues");
default
disp("Another solution status:")
disk(solsta)
end
else
warning("Solution does not exists");
end
end
