15 List of examples¶
List of examples shipped in the distribution of Rmosek package:
File |
Description |
|---|---|
A simple problem using affine conic constraints |
|
A simple problem using affine conic constraints |
|
A simple conic exponential problem |
|
A simple conic quadratic problem |
|
A simple geometric program (GP) in conic form |
|
A Hello World example |
|
A simple linear problem |
|
Implements logistic regression and simple log-sum-exp (CEO) |
|
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 |
|
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 quadratic problem |
|
Demonstrate how to modify and re-optimize a linear problem |
|
Demonstrates proper response handling |
|
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. |
|
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.
affco1.R
##
# File : affco1.R
#
# 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
#
##
library("Rmosek")
affco1 <- function()
{
prob <- list(sense="max")
# Variables [x1; x2; t1; t2]
prob$c <- c(0, 0, 1, 1)
# Linear inequality x_1 - x_2 <= 1
prob$A <- Matrix(c(1, -1, 0, 0), nrow=1, sparse=TRUE)
prob$bc <- rbind(blc=-Inf, buc=1)
prob$bx <- rbind(blx=rep(-Inf,4), bux=rep(Inf,4))
# The quadratic cone
FQ <- rbind(c(0,0,0,0), c(1,0,0,0), c(0,1,0,0))
gQ <- c(1, -0.5, -0.6)
cQ <- matrix(list("QUAD", 3, NULL), nrow=3, ncol=1)
# The power cone for (x_1, 1, t_1) \in POW3^(1/3,2/3)
FP1 <- rbind(c(1,0,0,0), c(0,0,0,0), c(0,0,1,0))
gP1 <- c(0, 1, 0)
cP1 <- matrix(list("PPOW", 3, c(1/3, 2/3)), nrow=3, ncol=1)
# The power cone for (x_1+x_2+0.1, 1, t_2) \in POW3^(1/4,3/4)
FP2 <- rbind(c(1,1,0,0), c(0,0,0,0), c(0,0,0,1))
gP2 <- c(0.1, 1, 0)
cP2 <- matrix(list("PPOW", 3, c(1.0, 3.0)), nrow=3, ncol=1)
# All cones
prob$F <- rbind(FQ, FP1, FP2)
prob$g <- cbind(gQ, gP1, gP2)
prob$cones <- cbind(cQ, cP1, cP2)
rownames(prob$cones) <- c("type","dim","conepar")
r <- mosek(prob, list(soldetail=1))
stopifnot(identical(r$response$code, 0))
print(r$sol$itr$pobjval)
print(r$sol$itr$xx[1:2])
}
affco1()
affco2.R
##
# File : affco2.R
#
# 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
#
##
library("Rmosek")
firstHittingTime <- function(n, z, a, d)
{
prob <- list(sense="min")
# Variables [t, y1, ..., yn]
prob$A <- Matrix(nrow=0, ncol=n+1)
prob$bx<- rbind(rep(-Inf,n+1), rep(Inf,n+1))
prob$c <- c(1, rep(0, n))
# Quadratic cone
FQ <- Diagonal(n+1, c(0, z))
gQ <- c(d, rep(0, n))
# All exponential cones
FE <- sparseMatrix( c(seq(1,3*n,by=3), seq(3,3*n,by=3)),
c(2:(n+1), rep(1,n)),
x=c(rep(1,n), t(a)))
gE <- rep(c(0, 1, 0), n)
# Assemble input data
prob$F <- rbind(FQ, FE)
prob$g <- c(gQ, gE)
prob$cones <- cbind(matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1),
matrix(list("PEXP", 3, NULL), nrow=3, ncol=n))
rownames(prob$cones) <- c("type","dim","conepar")
# Solve
r <- mosek(prob, list(soldetail=1))
stopifnot(identical(r$response$code, 0))
r$sol$itr$xx[1]
}
n = 3
z = c(3.3, 2.2, 1.3)
a = c(-0.08, -0.3, -0.06)
d = 0.5
t <- firstHittingTime(n, z, a, d)
print(sprintf("t = %.4e", t))
ceo1.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : ceo1.R
#
# Purpose : To demonstrate how to solve a small conic exponential
# optimization problem using the Rmosek package.
##
library("Rmosek")
ceo1 <- function()
{
# Specify the non-conic part of the problem.
prob <- list(sense="min")
prob$c <- c(1, 1, 0)
prob$A <- Matrix(c(1, 1, 1), nrow=1, sparse=TRUE)
prob$bc <- rbind(blc=1,
buc=1)
prob$bx <- rbind(blx=rep(-Inf,3),
bux=rep( Inf,3))
# Specify the affine conic constraints.
prob$F <- rbind(c(1,0,0),c(0,1,0),c(0,0,1))
prob$g <- c(0, 0, 0)
prob$cones <- matrix(list("PEXP", 3, NULL), nrow=3, ncol=1)
rownames(prob$cones) <- c("type","dim", "conepar")
# Solve the problem
r <- mosek(prob)
# Return the solution
stopifnot(identical(r$response$code, 0))
r$sol
}
ceo1()
cqo1.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : cqo1.R
#
# Purpose : To demonstrate how to solve a small conic quadratic
# optimization problem using the Rmosek package.
##
library("Rmosek")
cqo1 <- function()
{
# Specify the non-conic part of the problem.
prob <- list(sense="min")
prob$c <- c(0, 0, 0, 1, 1, 1)
prob$A <- Matrix(c(1, 1, 2, 0, 0, 0), nrow=1, sparse=TRUE)
prob$bc <- rbind(blc=1,
buc=1)
prob$bx <- rbind(blx=c(rep(0,3), rep(-Inf,3)),
bux=rep(Inf,6))
# Specify the affine conic constraints.
# NOTE: The F matrix is internally stored in the sparse
# triplet form. Use 'giveCsparse' or 'repr' option
# in the sparseMatrix() call to construct the F
# matrix directly in the sparse triplet form.
prob$F <- sparseMatrix(i=c(1, 2, 3, 4, 5, 6),
j=c(4, 1, 2, 5, 6, 3),
x=c(1, 1, 1, 1, 1, 1),
dims = c(6,6))
prob$g <- c(1:6)*0
prob$cones <- matrix(list(), nrow=3, ncol=2)
rownames(prob$cones) <- c("type","dim","conepar")
prob$cones[,1] <- list("QUAD", 3, NULL)
prob$cones[,2] <- list("RQUAD",3, NULL)
#
# Use cbind to extend this chunk of ACCs, if needed:
#
# oldcones <- prob$cones
# prob$cones <- cbind(oldcones, newcones)
#
# Solve the problem
r <- mosek(prob)
# Return the solution
stopifnot(identical(r$response$code, 0))
r$sol
}
cqo1()
gp1.R
#
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : gp1.R
#
# 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),
library("Rmosek")
gp1 <- function()
{
# Input data
Awall <- 200.0
Afloor <- 50.0
alpha <- 2.0
beta <- 10
gamma <- 2
delta <- 10
# Objective
prob <- list(sense="max")
prob$c <- 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 <- rbind(c(0, 0, 0, 1, 1),
c(0, 1, 1, 0, 0),
c(1,-1, 0, 0, 0),
c(0,-1, 1, 0, 0))
prob$bc <- rbind(c(1, -Inf, log(alpha), log(gamma)),
c(1, log(Afloor), log(beta), log(delta)))
prob$bx <- rbind(c(-Inf, -Inf, -Inf, -Inf, -Inf),
c( Inf, Inf, Inf, Inf, Inf))
# The conic part FX+g \in Kexp x Kexp
# x y z u v
# [ 0 0 0 1 0 ] 0
# [ 0 0 0 0 0 ] 1 in Kexp
# [ 1 1 0 0 0 ] log(2/Awall)
#
# [ 0 0 0 0 1 ] 0
# [ 0 0 0 0 0 ] 1 in Kexp
# [ 1 0 1 0 0 ] + log(2/Awall)
#
# NOTE: The F matrix is internally stored in the sparse
# triplet form. Use 'giveCsparse' or 'repr' option
# in the sparseMatrix() call to construct the F
# matrix directly in the sparse triplet form.
prob$F <- sparseMatrix(i = c(1, 3, 3, 4, 6, 6),
j = c(4, 1, 2, 5, 1, 3),
x = c(1, 1, 1, 1, 1, 1),
dims = c(6,5))
prob$g <- c(0, 1, log(2/Awall), 0, 1, log(2/Awall))
prob$cones <- matrix(rep(list("PEXP", 3, NULL), 2), ncol=2)
rownames(prob$cones) <- c("type","dim","conepar")
# Optimize and print results
r <- mosek(prob, list(soldetail=1))
print(exp(r$sol$itr$xx[1:3]))
}
gp1()
helloworld.R
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: helloworld.R
#
# The most basic example of how to get started with MOSEK.
library("Rmosek")
prob <- list(sense="min") # Minimization problem
prob$A <- Matrix(nrow=0, ncol=1) # 0 constraints, 1 variable
prob$bx <- rbind(blx=2.0, bux=3.0) # Bounds on the only variable
prob$c <- c(1.0) # The objective coefficient
# Optimize
r <- mosek(prob)
# Print answer
r$sol$itr$xx
lo1.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : lo1.R
#
# Purpose : To demonstrate how to solve a small linear
# optimization problem using the Rmosek package.
##
library("Rmosek")
lo1 <- function()
{
prob <- list()
# Objective sense (maximize or minimize)
prob$sense <- "max"
# Objective coefficients
prob$c <- c(3, 1, 5, 1)
# Specify matrix 'A' in sparse format.
asubi <- c(1, 1, 1, 2, 2, 2, 2, 3, 3)
asubj <- c(1, 2, 3, 1, 2, 3, 4, 2, 4)
aval <- c(3, 1, 2, 2, 1, 3, 1, 2, 3)
prob$A <- sparseMatrix(asubi,asubj,x=aval)
# Bound values for constraints
prob$bc <- rbind(blc=c(30, 15, -Inf),
buc=c(30, Inf, 25))
# Bound values for variables
prob$bx <- rbind(blx=rep(0,4),
bux=c(Inf, 10, Inf, Inf))
# Solve the problem
r <- mosek(prob)
# Return the solution
stopifnot(identical(r$response$code, 0))
r$sol
}
lo1()
logistic.R
##
# File : logistic.R
#
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# Description : Implements logistic regression with regulatization.
#
# Demonstrates using the exponential cone and log-sum-exp.
#
# Plots an example for 2D datasets.
#
##
library("Rmosek")
logisticRegression <- function(X, y, lamb)
{
prob <- list(sense="min")
n <- dim(X)[1];
d <- dim(X)[2];
# Variables: r, theta(d), t(n), z1(n), z2(n)
prob$c <- c(lamb, rep(0,d), rep(1, n), rep(0,n), rep(0,n));
prob$bx <-rbind(rep(-Inf,1+d+3*n), rep(Inf,1+d+3*n));
# z1 + z2 <= 1
prob$A <- sparseMatrix( rep(1:n, 2),
c((1:n)+1+d+n, (1:n)+1+d+2*n),
x = rep(1, 2*n));
prob$bc <- rbind(rep(-Inf, n), rep(1, n));
# (r, theta) \in \Q
FQ <- cbind(diag(rep(1, d+1)), matrix(0, d+1, 3*n));
gQ <- rep(0, 1+d);
# (z1(i), 1, -t(i)) \in \EXP,
# (z2(i), 1, (1-2y(i))*X(i,) - t(i)) \in \EXP
FE <- Matrix(nrow=0, ncol = 1+d+3*n);
for(i in 1:n) {
FE <- rbind(FE,
sparseMatrix( c(1, 3, 4, rep(6, d), 6),
c(1+d+n+i, 1+d+i, 1+d+2*n+i, 2:(d+1), 1+d+i),
x = c(1, -1, 1, (1-2*y[i])*X[i,], -1),
dims = c(6, 1+d+3*n) ) );
}
gE <- rep(c(0, 1, 0, 0, 1, 0), n);
prob$F <- rbind(FQ, FE)
prob$g <- c(gQ, gE)
prob$cones <- cbind(matrix(list("QUAD", 1+d, NULL), nrow=3, ncol=1),
matrix(list("PEXP", 3, NULL), nrow=3, ncol=2*n));
rownames(prob$cones) <- c("type","dim","conepar")
# Solve, no error handling!
r <- mosek(prob, list(soldetail=1))
# Return theta
r$sol$itr$xx[2:(d+1)]
}
mapFeature <- function(x, d)
{
g = as.matrix(expand.grid(0:d,0:d));
g = g[g[,1]+g[,2]<=d,];
apply(g, 1, function(gel) x[1]^gel[1]*x[2]^gel[2]);
}
mapFeatures <- function(X, d)
{
t(apply(X, 1, function(x) mapFeature(x, d)));
}
runExample <-function(x1, x2, y, deg, lamb, pngname)
{
X = mapFeatures(cbind(x1,x2),deg);
theta = logisticRegression(X, y, lamb);
print(theta);
xplot <- seq(-1, 1, 0.05);
yplot <- seq(-1, 1, 0.05);
nplot <- length(xplot);
zplot <- matrix(mapFeatures(as.matrix(expand.grid(xplot, yplot)),deg) %*% theta, ncol=nplot, nrow=nplot);
# Uncomment to plot
#png(pngname);
#plot(x1, x2, pch=19, col=y+1);
#contour(xplot, xplot, zplot, levels=c(0), add=TRUE);
}
example1 <- function()
{
# Example from documentation is contained in
# https://datascienceplus.com/wp-content/uploads/2017/02/ex2data2.txt
ex <- read.csv("ex2data2.txt");
x1 <- ex[,1];
x2 <- ex[,2];
y <- ex[,3];
deg <- 6;
for(lamb in c(1e-0,1e-2,1e-4))
{
runExample(x1, x2, y, deg, lamb, sprintf("ex2-%f.png", lamb));
}
}
example2 <- function()
{
# Example 2: discover a circle
n <- 300;
deg <- 2;
x1 <- runif(n, -1, 1);
x2 <- runif(n, -1, 1);
y <- ifelse(x1^2 + x2^2 - 0.69 < 0, 0, 1);
x1 <- x1 + runif(n, -0.2, 0.2);
x2 <- x2 + runif(n, -0.2, 0.2);
runExample(x1, x2, y, deg, 0.1, "circle.png");
}
#example1();
example2();
mico1.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : mico1.R
#
# Purpose : To demonstrate how to solve a small mixed integer conic
# optimization problem using the Rmosek package.
#
# minimize x^2 + y^2
# subject to x >= e^y + 3.8
# x, y - integer
##
library("Rmosek")
mico1 <- function()
{
# Specify the continuous part of the problem.
# Variables are [t; x; y]
prob <- list(sense="min")
prob$c <- c(1, 0, 0)
prob$A <- Matrix(nrow=0, ncol=3) # 0 constraints, 3 variables
prob$bx <- rbind(blx=rep(-Inf,3), bux=rep(Inf,3))
# Conic part of the problem
prob$F <- rbind(diag(3), c(0,1,0), c(0,0,0), c(0,0,1))
prob$g <- c(0, 0, 0, -3.8, 1, 0)
prob$cones <- cbind(matrix(list("QUAD", 3, NULL), nrow=3, ncol=1),
matrix(list("PEXP", 3, NULL), nrow=3, ncol=1))
rownames(prob$cones) <- c("type","dim","conepar")
# Specify the integer constraints
prob$intsub <- c(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 <- c(100, 9, -1)
prob$iparam <- list(MIO_CONSTRUCT_SOL=1)
# Solve the problem
r <- mosek(prob, list(getinfo=TRUE))
# Return the solution
stopifnot(identical(r$response$code, 0))
print(r$sol$int$xx[2:3])
}
mico1()
milo1.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : milo1.R
#
# Purpose : To demonstrate how to solve a small mixed integer linear
# optimization problem using the Rmosek package.
##
library("Rmosek")
milo1 <- function()
{
# Specify the continuous part of the problem.
prob <- list(sense="max")
prob$c <- c(1, 0.64)
prob$A <- Matrix(rbind(c(50, 31),
c( 3, -2)), sparse=TRUE)
prob$bc <- rbind(blc=c(-Inf, -4),
buc=c( 250, Inf))
prob$bx <- rbind(blx=c( 0, 0),
bux=c(Inf, Inf))
# Specify the integer constraints
prob$intsub <- c(1 ,2)
# Solve the problem
r <- mosek(prob)
# Return the solution
stopifnot(identical(r$response$code, 0))
r$sol
}
milo1()
mioinitsol.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : mioinitsol.R
#
# Purpose : To demonstrate how to solve a linear mixed-integer
# problem with a start guess.
##
library("Rmosek")
mioinitsol <- function()
{
# Specify the problem.
prob <- list(sense="max")
prob$c <- c(7, 10, 1, 5)
prob$A <- Matrix(c(1, 1, 1, 1), nrow=1, sparse=TRUE)
prob$bc <- rbind(blc=-Inf,
buc=2.5)
prob$bx <- rbind(blx=rep(0,4),
bux=rep(Inf,4))
prob$intsub <- c(1 ,2, 3)
# Specify start guess for the integer variables.
prob$sol$int$xx <- c(1, 1, 0, NaN)
# Request constructing the solution from integer variable values
prob$iparam <- list(MIO_CONSTRUCT_SOL=1)
# Solve the problem
r <- mosek(prob, list(getinfo=TRUE))
# The solution
stopifnot(identical(r$response$code, 0))
print(r$sol$int$xx)
# Was the initial solution used ?
print(r$iinfo$MIO_CONSTRUCT_SOLUTION)
print(r$dinfo$MIO_CONSTRUCT_SOLUTION_OBJ)
}
mioinitsol()
normex.R
##
#
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : normex.R
#
# Purpose: Demonstrates various norm minimization problems
#
# * least squares regression
# * ridge regularization
# * lasso regularization
# * p-norm minimization
#
##
library("Rmosek")
# This function runs a few examples
normex <- function()
{
# Create some random data
# We will attempt to fit various polynomials to
# a noisy degree 3 polynomial p(t)
N <- 100
x1 <- sort(runif(N))
p <- function(t) { 1.5*t^3-2*t^2+0.5*t^1-1 }
y1 <- p(x1) + 0.01*runif(N);
# Least squares regression
deg<- 5
xp <- outer(x1, 0:deg, '^')
x <- norm_lse(xp, y1, matrix(0, 0, deg+1), numeric(0))
print(sprintf("LSE: best fit %.2f x^%d + %.2f x^%d + ... + %.2f", x[deg+1], deg, x[deg], deg-1, x[1]))
# With ridge regularization
x <- norm_lse_reg(xp, y1, matrix(0, 0, deg+1), numeric(0), 0.2)
print(sprintf("Ridge: best fit %.2f x^%d + %.2f x^%d + ... + %.2f", x[deg+1], deg, x[deg], deg-1, x[1]))
# With squared version of ridge regularization
x <- norm_lse_reg_quad(xp, y1, matrix(0, 0, deg+1), numeric(0), 0.2)
print(sprintf("QRidge: best fit %.2f x^%d + %.2f x^%d + ... + %.2f", x[deg+1], deg, x[deg], deg-1, x[1]))
# Completely random large data for lasso example
N <- 10
K <- 3000
F <- matrix(rnorm(K*N), K, N)
g <- runif(K)
print(sprintf("Lasso regularization"))
for (gamma in c(0.01, 0.1, 0.3, 0.6, 0.9, 1.3))
{
x <- norm_lse_lasso(F, g, matrix(0,0,N), numeric(0), gamma)
z <- data.frame(x)
print(sprintf('Gamma %.4f density %.0f |Fx-g|_2: %.4f', gamma, sum(abs(z)>1e-6)/N*100, norm(F%*%x-g, type="2")))
}
# Example with the p-norm cone for various p
# We add a far outlier to the first example
x12 <- c(x1, 0.73)
y12 <- c(y1, -0.99)
xp2 <- outer(x12, 0:deg, '^')
for (p in c(1.1, 2.0, 3.0, 6.0))
{
x <- norm_p_norm(xp2, y12, matrix(0, 0, deg+1), numeric(0), p)
print(sprintf("p = %.1f: best fit %.2f x^%d + %.2f x^%d + ... + %.2f", p, x[deg+1], deg, x[deg], deg-1, x[1]))
}
}
# Least squares regression
# minimize \|Fx-g\|_2
norm_lse <- function(F,g,A,b)
{
n <- dim(F)[2]
k <- length(g)
m <- dim(A)[1]
# Linear constraints in [x; t]
prob <- list(sense="min")
prob$A <- cbind(A, rep(0, m))
prob$bx <- rbind(rep(-Inf, n+1), rep(Inf, n+1))
prob$bc <- rbind(b, b)
prob$c <- c(rep(0, n), 1)
# Affine conic constraint
prob$F <- rbind( c(rep(0,n), 1),
cbind(F, rep(0, k)) )
prob$g <- c(0, -g)
prob$cones <- matrix(list("QUAD", k+1, NULL))
rownames(prob$cones) <- c("type","dim","conepar")
# Solve
r <- mosek(prob, list(verbose=1))
stopifnot(identical(r$response$code, 0))
r$sol$itr$xx[1:n]
}
# Least squares regression with regularization
# minimize \|Fx-g\|_2 + gamma*\|x\|_2
norm_lse_reg <- function(F,g,A,b,gamma)
{
n <- dim(F)[2]
k <- length(g)
m <- dim(A)[1]
# Linear constraints in [x; t1; t2]
prob <- list(sense="min")
prob$A <- cbind(A, rep(0, m), rep(0, m))
prob$bx <- rbind(rep(-Inf, n+2), rep(Inf, n+2))
prob$bc <- rbind(b, b)
prob$c <- c(rep(0, n), 1, gamma)
# Affine conic constraint
prob$F <- rbind( c(rep(0,n), 1, 0),
cbind(F, rep(0, k), rep(0, k)),
c(rep(0,n), 0, 1),
cbind(diag(n), rep(0, n), rep(0, n)))
prob$g <- c(0, -g, rep(0, n+1))
prob$cones <- cbind(matrix(list("QUAD", k+1, NULL)),
matrix(list("QUAD", n+1, NULL)))
rownames(prob$cones) <- c("type","dim","conepar")
# Solve
r <- mosek(prob, list(verbose=1))
stopifnot(identical(r$response$code, 0))
r$sol$itr$xx[1:n]
}
# Least squares regression with regularization
# The "classical" quadratic version
# minimize \|Fx-g\|_2^2 + gamma*\|x\|_2^2
norm_lse_reg_quad <- function(F,g,A,b,gamma)
{
n <- dim(F)[2]
k <- length(g)
m <- dim(A)[1]
# Linear constraints in [x; t1; t2]
prob <- list(sense="min")
prob$A <- cbind(A, rep(0, m), rep(0, m))
prob$bx <- rbind(rep(-Inf, n+2), rep(Inf, n+2))
prob$bc <- rbind(b, b)
prob$c <- c(rep(0, n), 1, gamma)
# Affine conic constraint
prob$F <- rbind( c(rep(0,n), 1, 0),
c(rep(0, n+2)),
cbind(F, rep(0, k), rep(0, k)),
c(rep(0,n), 0, 1),
c(rep(0, n+2)),
cbind(diag(n), rep(0, n), rep(0, n)))
prob$g <- c(0, 0.5, -g, 0, 0.5, rep(0, n))
prob$cones <- cbind(matrix(list("RQUAD", k+2, NULL)),
matrix(list("RQUAD", n+2, NULL)))
rownames(prob$cones) <- c("type","dim","conepar")
# Solve
r <- mosek(prob, list(verbose=1))
stopifnot(identical(r$response$code, 0))
r$sol$itr$xx[1:n]
}
# Least squares regression with lasso regularization
# minimize \|Fx-g\|_2 + gamma*\|x\|_1
norm_lse_lasso <- function(F,g,A,b,gamma)
{
n <- dim(F)[2]
k <- length(g)
m <- dim(A)[1]
# Linear constraints in [x; u; t1; t2]
prob <- list(sense="min")
prob$A <- rbind(cbind(A, matrix(0, m, n+2)),
cbind(diag(n), diag(n), matrix(0, n, 2)),
cbind(-diag(n), diag(n), matrix(0, n, 2)),
c(rep(0, n), rep(-1, n), 0, 1))
prob$bx <- rbind(rep(-Inf, 2*n+2), rep(Inf, 2*n+2))
prob$bc <- rbind(c(b, rep(0, 2*n+1)), c(b, rep(Inf, 2*n+1)))
prob$c <- c(rep(0, 2*n), 1, gamma)
# Affine conic constraint
prob$F <- rbind(c(rep(0, 2*n), 1, 0),
cbind(F, matrix(0, k, n+2)))
prob$g <- c(0, -g)
prob$cones <- matrix(list("QUAD", k+1, NULL))
rownames(prob$cones) <- c("type","dim","conepar")
# Solve
r <- mosek(prob, list(verbose=1))
stopifnot(identical(r$response$code, 0))
r$sol$itr$xx[1:n]
}
# P-norm minimization
# minimize \|Fx-g\|_p
norm_p_norm <- function(F,g,A,b,p)
{
n <- dim(F)[2]
k <- length(g)
m <- dim(A)[1]
# Linear constraints in [x; r; t]
prob <- list(sense="min")
prob$A <- rbind(cbind(A, matrix(0, m, k+1)),
c(rep(0, n), rep(1, k), -1))
prob$bx <- rbind(rep(-Inf, n+k+1), rep(Inf, n+k+1))
prob$bc <- rbind(c(b, 0), c(b, 0))
prob$c <- c(rep(0, n+k), 1)
# Permutation matrix which picks triples (r_i, t, F_ix-g_i)
M <- cbind(sparseMatrix(seq(1,3*k,3), seq(1,k), x=rep(1,k), dims=c(3*k,k)),
sparseMatrix(seq(2,3*k,3), seq(1,k), x=rep(1,k), dims=c(3*k,k)),
sparseMatrix(seq(3,3*k,3), seq(1,k), x=rep(1,k), dims=c(3*k,k)))
# Affine conic constraint
prob$F <- M %*% rbind(cbind(matrix(0, k, n), diag(k), rep(0, k)),
cbind(matrix(0, k, n+k), rep(1, k)),
cbind(F, matrix(0, k, k+1)))
prob$g <- as.numeric(M %*% c(rep(0, 2*k), -g))
prob$cones <- matrix(list("PPOW", 3, c(1, p-1)), nrow=3, ncol=k)
rownames(prob$cones) <- c("type","dim","conepar")
# Solve
r <- mosek(prob, list(verbose=1))
stopifnot(identical(r$sol$itr$prosta, "PRIMAL_AND_DUAL_FEASIBLE"))
r$sol$itr$xx[1:n]
}
normex()
parameters.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : parameters.R
#
# Purpose : To demonstrate how to set parameters
# using the Rmosek package.
##
library("Rmosek")
parameters <- function()
{
# Specify the optimization problem
prob <- list(sense="max")
prob$c <- c(3, 1, 5, 1)
prob$A <- sparseMatrix(i=c(1, 1, 1, 2, 2, 2, 2, 3, 3),
j=c(1, 2, 3, 1, 2, 3, 4, 2, 4),
x=c(3, 1, 2, 2, 1, 3, 1, 2, 3))
prob$bc <- rbind(blc=c(30, 15, -Inf),
buc=c(30, Inf, 25))
prob$bx <- rbind(blx=rep(0,4),
bux=c(Inf, 10, Inf, Inf))
# Specify the parameters
# Set log level (integer parameter)
# and select interior-point optimizer... (integer parameter)
# ... without basis identification (integer parameter)
prob$iparam <- list(LOG=0, OPTIMIZER="OPTIMIZER_INTPNT", INTPNT_BASIS="BI_NEVER")
# Set relative gap tolerance (double parameter)
prob$dparam <- list(INTPNT_CO_TOL_REL_GAP=1.0e-7)
# Solve the problem
r <- mosek(prob)
# Demonstrate retrieving information items
opts <- list(getinfo=TRUE)
r <- mosek(prob, opts)
print(r$dinfo$OPTIMIZER_TIME)
print(r$iinfo$INTPNT_ITER)
# Return the solution
stopifnot(identical(r$response$code, 0))
r$sol
}
parameters()
pinfeas.R
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: pinfeas.R
#
# Purpose: Demonstrates how to fetch a primal infeasibility certificate
# for a linear problem
#
library("Rmosek")
# Set up a simple linear problem from the manual for test purposes
testProblem <- function()
{
prob <- list(sense="min")
prob$c <- c(1, 2, 5, 2, 1, 2, 1)
prob$A <- sparseMatrix(i=c(1,1,2,2,3,3,3,4,4,5,6,6,7,7),
j=c(1,2,3,4,5,6,7,1,5,2,3,6,4,7),
x=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1),
dims = c(7,7))
prob$bc <- rbind(blc=c(-Inf, -Inf, -Inf, 1100, 200, 500, 500),
buc=c(200, 1000, 1000, 1100, 200, 500, 500))
prob$bx <- rbind(blx=c(0, 0, 0, 0, 0, 0, 0),
bux=c(Inf, Inf, Inf, Inf, Inf , Inf ,Inf))
prob
}
# Analyzes and prints infeasibility contributing elements
analyzeCertificate <- function(sl, # dual values for lower bounds
su, # dual values for upper bounds
eps) # tolerance determining when a dual value is considered important
{
n <- length(sl)
for(i in 1:n) {
if (abs(sl[i]) > eps) print(sprintf("#%d: lower, dual = %e", i, sl[i]))
if (abs(su[i]) > eps) print(sprintf("#%d: upper, dual = %e", i, su[i]))
}
}
pinfeas <- function()
{
# In this example we set up a simple problem
prob <- testProblem()
# Perform the optimization.
r <- mosek(prob)
# Use the line below instead to disable log output
#r <- mosek(prob, list(verbose=0))
# Check problem status
if (r$sol$itr$prosta == 'PRIMAL_INFEASIBLE') {
# Set the tolerance at which we consider a dual value as essential
eps <- 1e-7
print("Variable bounds important for infeasibility: ")
analyzeCertificate(r$sol$itr$slx, r$sol$itr$sux, eps)
print("Constraint bounds important for infeasibility: ")
analyzeCertificate(r$sol$itr$slc, r$sol$itr$suc, eps)
}
else {
print("The problem is not primal infeasible, no certificate to show")
}
}
pinfeas()
portfolio_1_basic.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_1_basic.R
#
# Purpose : To implement a basic Markowitz optimization model computing
# the optimal expected return and portfolio for a given risk.
##
library("Rmosek")
BasicMarkowitz <- function(
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
{
prob <- list(sense="max")
prob$c <- mu
prob$A <- Matrix(1.0, ncol=n)
prob$bc <- rbind(blc=w+sum(x0),
buc=w+sum(x0))
prob$bx <- rbind(blx=rep(0.0,n),
bux=rep(Inf,n))
# Specify the affine conic constraints.
NUMCONES <- 1
prob$F <- rbind(
Matrix(0.0,ncol=n),
GT
)
prob$g <- c(gamma,rep(0,n))
prob$cones <- matrix(list(), nrow=3, ncol=NUMCONES)
rownames(prob$cones) <- c("type","dim","conepar")
prob$cones[-3,1] <- list("QUAD", n+1)
# Solve the problem
r <- mosek(prob,list(verbose=1))
stopifnot(identical(r$response$code, 0))
# Check if the interior point solution is an optimal point
# See https://docs.mosek.com/latest/rmosek/accessing-solution.html about handling solution statuses.
stopifnot(identical(r$sol$itr$solsta, 'OPTIMAL'))
# Return the solution
x <- r$sol$itr$xx
list(expret=drop(mu %*% x), stddev=gamma, x=x)
}
# Example of input
n <- 8
w <- 59.0
mu <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0 <- c(8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0)
gamma <- 36.0
GT <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
c(0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
c(0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
c(0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
c(0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ),
c(0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187),
c(0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327),
c(0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ) )
r <- BasicMarkowitz(n,mu,GT,x0,w,gamma)
print(sprintf('Expected return: %.4e Std. deviation: %.4e', r$expret, r$stddev))
portfolio_2_frontier.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_2_frontier.R
#
# Purpose : To implements a basic portfolio optimization model.
# Computes points on the efficient frontier.
##
library("Rmosek")
EfficientFrontier <- function(
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 risk penalties (we maximize expected return - alpha * variance)
{
prob <- list(sense="max")
prob$A <- cbind(Matrix(1.0, ncol=n), 0.0)
prob$bc <- rbind(blc=w+sum(x0),
buc=w+sum(x0))
prob$bx <- rbind(blx=c(rep(0.0,n), -Inf),
bux=rep(Inf,n+1))
# Specify the affine conic constraints.
NUMCONES <- 1
prob$F <- rbind(
cbind(Matrix(0.0,ncol=n), 1.0),
rep(0, n+1),
cbind(GT , 0.0)
)
prob$g <- c(0, 0.5, rep(0, n))
prob$cones <- matrix(list(), nrow=3, ncol=NUMCONES)
rownames(prob$cones) <- c("type","dim","conepar")
prob$cones[-3,1] <- list("RQUAD", n+2)
frontier <- matrix(NaN, ncol=3, nrow=length(alphas))
colnames(frontier) <- c("alpha","exp.ret.","variance")
for (i in seq_along(alphas))
{
prob$c <- c(mu, -alphas[i])
r <- mosek(prob, list(verbose=1))
stopifnot(identical(r$response$code, 0))
# Check if the interior point solution is an optimal point
# See https://docs.mosek.com/latest/rmosek/accessing-solution.html about handling solution statuses.
stopifnot(identical(r$sol$itr$solsta, 'OPTIMAL'))
x <- r$sol$itr$xx[1:n]
gamma <- r$sol$itr$xx[n+1]
frontier[i,] <- c(alphas[i], drop(mu%*%x), gamma)
}
frontier
}
# Example of input
n <- 8
w <- 1.0
mu <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0 <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
GT <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
c(0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
c(0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
c(0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
c(0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ),
c(0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187),
c(0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327),
c(0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ) )
alphas <- c(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)
frontier <- EfficientFrontier(n,mu,GT,x0,w,alphas)
print(frontier)
portfolio_3_impact.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_3_impact.R
#
# Purpose : Implements a basic portfolio optimization model
# with x^(3/2) market impact costs.
##
library("Rmosek")
MarkowitzWithMarketImpact <- function(
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) # Market impacts (we use m_j|x_j-x0_j|^3/2 for j'th asset)
{
prob <- list(sense="max")
prob$c <- c(mu, rep(0,n))
prob$A <- cbind(Matrix(1.0,ncol=n), t(m))
prob$bc <- rbind(blc=w+sum(x0),
buc=w+sum(x0))
prob$bx <- rbind(blx=rep(0.0,2*n),
bux=rep(Inf,2*n))
# Specify the affine conic constraints.
# 1) Risk
Fr <- rbind(
Matrix(0.0,nrow=1,ncol=2*n),
cbind(GT, Matrix(0.0,nrow=n,ncol=n))
)
gr <- c(gamma,rep(0,n))
Kr <- matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1)
# 2) Market impact (t_j >= |x_j-x0_j|^3/2)
# [ t_j ]
# [ 1 ] \in PPOW(2,1)
# [ x_j - x0_j ]
Fm <- sparseMatrix(
i=c(seq(from=1,by=3,len=n), seq(from=3,by=3,len=n)),
j=c(seq(from=n+1,len=n), seq(from=1,len=n)),
x=c(rep(1.0,n), rep(1.0,n)),
dims=c(3*n, 2*n))
gm <- rep(c(0,1,0), n)
gm[seq(from=3,by=3,len=n)] <- -x0
Km <- matrix(list("PPOW", 3, c(2,1)), nrow=3, ncol=n)
prob$F <- rbind(Fr, Fm)
prob$g <- c(gr, gm)
prob$cones <- cbind(Kr, Km)
rownames(prob$cones) <- c("type","dim","conepar")
# Solve the problem
r <- mosek(prob,list(verbose=1))
stopifnot(identical(r$response$code, 0))
# Check if the interior point solution is an optimal point
# See https://docs.mosek.com/latest/rmosek/accessing-solution.html about handling solution statuses.
stopifnot(identical(r$sol$itr$solsta, 'OPTIMAL'))
# Return the solution
x <- r$sol$itr$xx[1:n]
tx <- r$sol$itr$xx[(n+1):(2*n)]
list(expret=drop(mu %*% x), stddev=gamma, cost=drop(m %*% tx), x=x)
}
# Example of input
n <- 8
w <- 1.0
mu <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0 <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
GT <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
c(0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
c(0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
c(0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
c(0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ),
c(0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187),
c(0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327),
c(0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ) )
gamma <- 0.36
m <- rep(0.01, n)
r <- MarkowitzWithMarketImpact(n,mu,GT,x0,w,gamma,m)
print(sprintf('Expected return: %.4e Std. deviation: %.4e Market impact cost: %.4e', r$expret, r$stddev, r$cost))
portfolio_4_transcost.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_4_transcost.R
#
# Purpose : Implements a basic portfolio optimization model
# with fixed setup costs and transaction costs
# as a mixed-integer problem.
##
library("Rmosek")
MarkowitzWithTransactionCosts <- function(
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, # Fixed transaction cost
g) # Linear part of transaction cost
{
# Upper bound on the traded amount
u <- w+sum(x0)
prob <- list(sense="max")
prob$c <- c(mu, rep(0,2*n))
# Specify 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 <- rbind(cbind(Matrix(1.0,ncol=n), t(g), t(f)),
cbind(Diagonal(n, 1.0), -Diagonal(n, 1.0), Matrix(0,n,n)),
cbind(Diagonal(n, 1.0), Diagonal(n, 1.0), Matrix(0,n,n)),
cbind(Matrix(0,n,n), Diagonal(n, 1.0), Diagonal(n, -u)))
prob$bc <- rbind(blc=c(w+sum(x0), rep(-Inf,n), x0, rep(-Inf,n)),
buc=c(w+sum(x0), x0, rep(Inf,n), rep(0.0,n)))
# No shortselling and the linear bound 0 <= y <= 1
prob$bx <- rbind(blx=c(rep(0.0,n), rep(-Inf,n), rep(0.0,n)),
bux=c(rep(Inf,n), rep(Inf, n), rep(1.0,n)))
# Specify the affine conic constraints for risk
prob$F <- rbind(
Matrix(0.0,nrow=1,ncol=3*n),
cbind(GT, Matrix(0.0,nrow=n,ncol=2*n))
)
prob$g <- c(gamma,rep(0,n))
prob$cones <- matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1)
rownames(prob$cones) <- c("type","dim","conepar")
# Demand y to be integer (hence binary)
prob$intsub <- (2*n+1):(3*n);
# Solve the problem
r <- mosek(prob,list(verbose=1))
stopifnot(identical(r$response$code, 0))
# Check if the integer solution is an optimal point
# See https://docs.mosek.com/latest/rmosek/accessing-solution.html about handling solution statuses.
stopifnot(identical(r$sol$int$solsta, 'INTEGER_OPTIMAL'))
# Return the solution
x <- r$sol$int$xx[1:n]
z <- r$sol$int$xx[(n+1):(2*n)]
y <- r$sol$int$xx[(2*n+1):(3*n)]
list(expret=drop(mu %*% x), stddev=gamma, cost = drop(f %*% y)+drop(g %*% z), x=x)
}
# Example of input
n <- 8
w <- 1.0
mu <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0 <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
GT <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
c(0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
c(0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
c(0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
c(0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ),
c(0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187),
c(0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327),
c(0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ) )
gamma <- 0.36
f <- rep(0.01,n)
g <- rep(0.001,n)
r <- MarkowitzWithTransactionCosts(n,mu,GT,x0,w,gamma,f,g)
print(sprintf("Markowitz portfolio optimization with transactions cost"))
print(sprintf("Expected return: %.4e Std. deviation: %.4e Transactions cost: %.4e", r$expret, r$stddev, r$cost))
portfolio_5_card.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_5_card.R
#
# Description : Implements a basic portfolio optimization model
# with cardinality constraints on number of assets traded.
##
library("Rmosek")
MarkowitzWithCardinality <- function(
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
k) # Cardinality bound
{
# Upper bound on the traded amount
u <- w+sum(x0)
prob <- list(sense="max")
prob$c <- c(mu, rep(0,2*n))
# Specify 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 <- rbind(cbind(Matrix(1.0,ncol=n), Matrix(0.0,ncol=2*n)),
cbind(Diagonal(n, 1.0), -Diagonal(n, 1.0), Matrix(0,n,n)),
cbind(Diagonal(n, 1.0), Diagonal(n, 1.0), Matrix(0,n,n)),
cbind(Matrix(0,n,n), Diagonal(n, 1.0), Diagonal(n, -u)),
cbind(Matrix(0.0,ncol=2*n), Matrix(1.0,ncol=n)))
prob$bc <- rbind(blc=c(w+sum(x0), rep(-Inf,n), x0, rep(-Inf,n), 0.0),
buc=c(w+sum(x0), x0, rep(Inf,n), rep(0.0,n), k))
# No shortselling and the linear bound 0 <= y <= 1
prob$bx <- rbind(blx=c(rep(0.0,n), rep(-Inf,n), rep(0.0,n)),
bux=c(rep(Inf,n), rep(Inf, n), rep(1.0,n)))
# Specify the affine conic constraints for risk
prob$F <- rbind(
Matrix(0.0,nrow=1,ncol=3*n),
cbind(GT, Matrix(0.0,nrow=n,ncol=2*n))
)
prob$g <- c(gamma,rep(0,n))
prob$cones <- matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1)
rownames(prob$cones) <- c("type","dim","conepar")
# Demand y to be integer (hence binary)
prob$intsub <- (2*n+1):(3*n);
# Solve the problem
r <- mosek(prob,list(verbose=1))
stopifnot(identical(r$response$code, 0))
# Check if the integer solution is an optimal point
# See https://docs.mosek.com/latest/rmosek/accessing-solution.html about handling solution statuses.
stopifnot(identical(r$sol$int$solsta, 'INTEGER_OPTIMAL'))
# Return the solution
x <- r$sol$int$xx[1:n]
list(card=k, expret=drop(mu %*% x), x=x)
}
# Example of input
n <- 8
w <- 1.0
mu <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0 <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
GT <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
c(0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
c(0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
c(0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
c(0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ),
c(0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187),
c(0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327),
c(0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ) )
gamma <- 0.25
print(sprintf("Markowitz portfolio optimization with cardinality constraints"))
for (k in 1:n) {
r <- MarkowitzWithCardinality(n,mu,GT,x0,w,gamma,k)
print(sprintf("Bound: %d Expected return: %.4e Solution:", r$card, r$expret))
print(r$x)
}
portfolio_6_factor.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_6_factor.R
#
# Purpose : Implements a basic Markowitz portfolio model of factor type
##
library("Rmosek")
BasicMarkowitz <- function(
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
{
k <- dim(G_factor_T)[1]
prob <- list(sense="max")
prob$c <- mu
prob$A <- Matrix(1.0, ncol=n)
prob$bc <- rbind(blc=w+sum(x0),
buc=w+sum(x0))
prob$bx <- rbind(blx=rep(0.0,n),
bux=rep(Inf,n))
# Specify the affine conic constraints.
NUMCONES <- 1
prob$F <- rbind(
Matrix(0.0,ncol=n),
G_factor_T,
diag(sqrt(theta))
)
prob$g <- c(gamma,rep(0,k+n))
prob$cones <- matrix(list(), nrow=3, ncol=NUMCONES)
rownames(prob$cones) <- c("type","dim","conepar")
prob$cones[-3,1] <- list("QUAD", k+n+1)
# Solve the problem
r <- mosek(prob,list(verbose=1))
stopifnot(identical(r$response$code, 0))
# Check if the interior point solution is an optimal point
# See https://docs.mosek.com/latest/rmosek/accessing-solution.html about handling solution statuses.
stopifnot(identical(r$sol$itr$solsta, 'OPTIMAL'))
# Return the solution
x <- r$sol$itr$xx
list(expret=drop(mu %*% x), stddev=gamma, x=x)
}
# Example of input
n <- 8
w <- 1.0
mu <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0 <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
gammas <- c(0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48)
B <- rbind(c(0.4256, 0.1869),
c(0.2413, 0.3877),
c(0.2235, 0.3697),
c(0.1503, 0.4612),
c(1.5325, -0.2633),
c(1.2741, -0.2613),
c(0.6939, 0.2372),
c(0.5425, 0.2116) )
S_F <- rbind(c(0.0620, 0.0577),
c(0.0577, 0.0908))
theta <- c(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 <- t(chol(S_F))
G_factor <- B %*% P
G_factor_T <- t(G_factor)
for (gamma in gammas)
{
r <- BasicMarkowitz(n,mu,G_factor_T,theta,x0,w,gamma)
print(sprintf('Gamma: %.4e Expected return: %.4e Std. deviation: %.4e', gamma, r$expret, r$stddev))
}
pow1.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : pow1.R
#
# Purpose : To demonstrate how to solve a small power cone
# optimization problem using the Rmosek package.
##
library("Rmosek")
pow1 <- function()
{
# Specify the non-conic part of the problem.
prob <- list(sense="max")
prob$c <- c(-1, 0, 0, 1, 1)
prob$A <- Matrix(c(1, 1, 0.5, 0, 0), nrow=1, sparse=TRUE)
prob$bc <- rbind(blc=2,
buc=2)
prob$bx <- rbind(blx=c(rep(-Inf,5)),
bux=c(rep( Inf,5)))
# Specify the affine conic constraints.
# NOTE: The F matrix is internally stored in the sparse
# triplet form. Use 'giveCsparse' or 'repr' option
# in the sparseMatrix() call to construct the F
# matrix directly in the sparse triplet form.
prob$F <- sparseMatrix(i=c(1, 2, 3, 4, 6),
j=c(1, 2, 4, 3, 5),
x=c(1, 1, 1, 1, 1),
dims = c(6,5))
prob$g <- c(0, 0, 0, 0, 1, 0)
prob$cones <- matrix(list(), nrow=3, ncol=2)
rownames(prob$cones) <- c("type","dim","conepar")
prob$cones[,1] <- list("PPOW", 3, c(0.2, 0.8))
prob$cones[,2] <- list("PPOW", 3, c(0.4, 0.6))
# Solve the problem
r <- mosek(prob)
# Return the solution
stopifnot(identical(r$response$code, 0))
r$sol
}
pow1()
qo1.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : qo1.R
#
# Purpose : To demonstrate how to solve a small quadratic
# optimization problem using the Rmosek package.
##
library("Rmosek")
qo1 <- function()
{
# Specify the non-quadratic part of the problem.
prob <- list(sense="min")
prob$c <- c(0, -1, 0)
prob$A <- Matrix(c(1, 1, 1), nrow=1, sparse=TRUE)
prob$bc <- rbind(blc=1,
buc=Inf)
prob$bx <- rbind(blx=rep(0,3),
bux=rep(Inf,3))
# Specify the quadratic objective matrix in triplet form.
prob$qobj$i <- c(1, 3, 2, 3)
prob$qobj$j <- c(1, 1, 2, 3)
prob$qobj$v <- c(2, -1, 0.2, 2)
# Solve the problem
r <- mosek(prob)
# Return the solution
stopifnot(identical(r$response$code, 0))
r$sol
}
qo1()
reoptimization.R
#
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : reoptimization.R
#
# Purpose : Demonstrates how to modify and re-optimize a linear problem
#
library("Rmosek")
reoptimization <- function()
{
# Specify the c vector.
prob <- list(sense="max")
prob$c <- c(1.5, 2.5, 3.0)
# Specify a in sparse format.
subi <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
subj <- c(1, 2, 3, 1, 2, 3, 1, 2, 3)
valij <- c(2, 4, 3, 3, 2, 3, 2, 3, 2)
prob$A <- sparseMatrix(subi,subj,x=valij);
# Specify bounds of the constraints.
prob$bc<- rbind(blc=rep(-Inf, 3),
buc=c(100000, 50000, 60000))
# Specify bounds of the variables.
prob$bx<- rbind(blx=rep(0,3),
bux=rep(Inf,3))
# Perform the optimization.
prob$iparam <- list(OPTIMIZER="OPTIMIZER_DUAL_SIMPLEX")
r <- mosek(prob)
# Show the optimal x solution.
print(r$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 <- list(bas=r$sol$bas)
r <- mosek(prob)
print(r$sol$bas$xx)
#####################################################
# Add a variable
prob$c <- c(prob$c, 1.0)
prob$A <- cbind(prob$A, c(4.0, 0.0, 1.0))
prob$bx <- cbind(prob$bx, c(0.0,Inf))
# Reoptimize with a new variable and hot-start
# All parts of the solution must be extended to the new dimensions.
prob$sol <- list(bas=r$sol$bas)
prob$sol$bas$xx <- c(prob$sol$bas$xx, 0.0)
prob$sol$bas$slx <- c(prob$sol$bas$slx, 0.0)
prob$sol$bas$sux <- c(prob$sol$bas$sux, 0.0)
prob$sol$bas$skx <- c(prob$sol$bas$skx, "UN")
r <- mosek(prob)
print(r$sol$bas$xx)
#####################################################
# Add a constraint
prob$A <- rbind(prob$A, c(1.0, 2.0, 1.0, 1.0))
prob$bc <- cbind(prob$bc, c(-Inf, 30000.0))
# Reoptimization with a new constraint
# All parts of the solution must be extended to the new dimensions.
prob$sol <- list(bas=r$sol$bas)
prob$sol$bas$xc <- c(prob$sol$bas$xc, 0.0)
prob$sol$bas$slc <- c(prob$sol$bas$slc, 0.0)
prob$sol$bas$suc <- c(prob$sol$bas$suc, 0.0)
prob$sol$bas$skc <- c(prob$sol$bas$skc, "UN")
r <- mosek(prob)
print(r$sol$bas$xx)
#####################################################
# Change constraint bounds
prob$bc<- rbind(blc=rep(-Inf, 4),
buc=c(80000, 40000, 50000, 22000))
r <- mosek(prob)
print(r$sol$bas$xx)
}
reoptimization()
response.R
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: response.R
#
# Purpose: This example demonstrates proper response handling
# for problems solved with the interior-point optimizers.
#
library("Rmosek")
# Set up a simple conic problem for test purposes
setupProblem <- function()
{
prob <- list(sense="min")
prob$c <- c(0, 0, 0, 1, 1, 1)
prob$A <- Matrix(c(1, 1, 2, 0, 0, 0), nrow=1, sparse=TRUE)
prob$bc <- rbind(blc=1, buc=1)
prob$bx <- rbind(blx=c(rep(0,3), rep(-Inf,3)), bux=rep(Inf,6))
prob$F <- sparseMatrix(i=c(1, 2, 3, 4, 5, 6),
j=c(4, 1, 2, 5, 6, 3),
x=c(1, 1, 1, 1, 1, 1),
dims = c(6,6))
prob$g <- c(1:6)*0
prob$cones <- matrix(list(), nrow=3, ncol=2)
rownames(prob$cones) <- c("type","dim","conepar")
prob$cones[,1] <- list("QUAD", 3, NULL)
prob$cones[,2] <- list("RQUAD", 3, NULL)
prob
}
response <- function()
{
# In this example we set up a simple problem
prob <- setupProblem()
# (Optionally) Uncomment the next line to get solution status Unknown
# prob$iparam <- list(INTPNT_MAX_ITERATIONS=1)
# Perform the optimization.
r <- mosek(prob, list(verbose=0))
# Use the line below instead to get more log output
#r <- mosek(prob, list(verbose=10))
# Expected result: The solution status of the interior-point solution is optimal.
# Check if there was a fatal error
if (r$response$code != 0 && startsWith(r$response$msg, "MSK_RES_ERR"))
{
print(sprintf("Optimization error: %s (%d),", r$response$msg, r$response$code))
}
else
{
if (r$sol$itr$solsta == 'OPTIMAL')
{
print('An optimal interior-point solution is located:')
print(r$sol$itr$xx)
}
else if (r$sol$itr$solsta == 'DUAL_INFEASIBLE_CER')
{
print('Dual infeasibility certificate found.')
}
else if (r$sol$itr$solsta == 'PRIMAL_INFEASIBLE_CER')
{
print('Primal infeasibility certificate found.')
}
else if (r$sol$itr$solsta == 'UNKNOWN')
{
# The solutions status is unknown. The termination code
# indicates why the optimizer terminated prematurely.
print('The solution status is unknown.')
print(sprintf('Termination code: %s (%d).', r$response$msg, r$response$code))
}
else
{
printf('An unexpected solution status is obtained.')
}
}
}
response()
sdo1.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : sdo1.R
#
# Purpose : To solve 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)
##
library("Rmosek")
getbarvarMatrix <- function(barvar, bardim)
{
N <- as.integer(bardim)
new("dspMatrix", x=barvar, uplo="L", Dim=c(N,N))
}
sdo1 <- function()
{
# Specify the non-matrix variable part of the problem.
prob <- list(sense="min")
prob$c <- c(1, 0, 0)
prob$A <- sparseMatrix(i=c(1, 2, 2),
j=c(1, 2, 3),
x=c(1, 1, 1), dims=c(2, 3))
prob$bc <- rbind(blc=c(1, 0.5),
buc=c(1, 0.5))
prob$bx <- rbind(blx=rep(-Inf,3),
bux=rep( Inf,3))
# NOTE: The F matrix is internally stored in the sparse
# triplet form. Use 'giveCsparse' or 'repr' option
# in the sparseMatrix() call to construct the F
# matrix directly in the sparse triplet form.
prob$F <- sparseMatrix(i=c(1,2,3),
j=c(1,2,3),
x=c(1,1,1),
dims = c(3,3))
prob$g <- c(1:3)*0
prob$cones <- cbind(list("QUAD", 3, NULL))
# Specify semidefinite matrix variables (one 3x3 block)
prob$bardim <- c(3)
# Block triplet format specifying the lower triangular part
# of the symmetric coefficient matrix 'barc':
prob$barc$j <- c(1, 1, 1, 1, 1)
prob$barc$k <- c(1, 2, 3, 2, 3)
prob$barc$l <- c(1, 2, 3, 1, 2)
prob$barc$v <- c(2, 2, 2, 1, 1)
# Block triplet format specifying the lower triangular part
# of the symmetric coefficient matrix 'barA':
prob$barA$i <- c(1, 1, 1, 2, 2, 2, 2, 2, 2)
prob$barA$j <- c(1, 1, 1, 1, 1, 1, 1, 1, 1)
prob$barA$k <- c(1, 2, 3, 1, 2, 3, 2, 3, 3)
prob$barA$l <- c(1, 2, 3, 1, 2, 3, 1, 1, 2)
prob$barA$v <- c(1, 1, 1, 1, 1, 1, 1, 1, 1)
# Solve the problem
r <- mosek(prob)
# Print matrix variable and return the solution
stopifnot(identical(r$response$code, 0))
print( list(barx=getbarvarMatrix(r$sol$itr$barx[[1]], prob$bardim[1])) )
r$sol
}
sdo1()
sdo2.R
#
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : sdo2.R
#
# 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.
library("Rmosek")
getbarvarMatrix <- function(barvar, bardim)
{
N <- as.integer(bardim)
new("dspMatrix", x=barvar, uplo="L", Dim=c(N,N))
}
sdo2 <- function()
{
# Sample data in sparse lower-triangular triplet form
C1_k <- c(1, 3);
C1_l <- c(1, 3);
C1_v <- c(1, 6);
A1_k <- c(1, 3, 3);
A1_l <- c(1, 1, 3);
A1_v <- c(1, 1, 2);
C2_k <- c(1, 2, 2, 3);
C2_l <- c(1, 1, 2, 3);
C2_v <- c(1, -3, 2, 1);
A2_k <- c(2, 2, 4);
A2_l <- c(1, 2, 4);
A2_v <- c(1, -1, -3);
b <- 23.0;
k <- -3.0;
# Specify the dimensions of the problem
prob <- list(sense="min");
# Two constraints
prob$A <- Matrix(nrow=2, ncol=0); # 2 constraints
prob$c <- numeric(0);
prob$bx <- rbind(blc=numeric(0),
buc=numeric(0));
# Dimensions of semidefinite matrix variables
prob$bardim <- c(3, 4);
# Constraint bounds
prob$bc <- rbind(blc=c(b, -Inf),
buc=c(b, k));
# Block triplet format specifying the lower triangular part
# of the symmetric coefficient matrix 'barc':
prob$barc$j <- c(rep(1, length(C1_v)),
rep(2, length(C2_v))); # Which PSD variable (j)
prob$barc$k <- c(C1_k, C2_k); # Entries: (k,l)->v
prob$barc$l <- c(C1_l, C2_l);
prob$barc$v <- c(C1_v, C2_v);
# Block triplet format specifying the lower triangular part
# of the symmetric coefficient matrix 'barA':
prob$barA$i <- c(rep(1, length(A1_v)+length(A2_v)), 2); # Which constraint (i)
prob$barA$j <- c(rep(1, length(A1_v)),
rep(2, length(A2_v)),
2); # Which PSD variable (j)
prob$barA$k <- c(A1_k, A2_k, 2); # Entries: (k,l)->v
prob$barA$l <- c(A1_l, A2_l, 1);
prob$barA$v <- c(A1_v, A2_v, 0.5);
# Solve the problem
r <- mosek(prob);
# Print matrix variable and return the solution
stopifnot(identical(r$response$code, 0));
print( list(X1=getbarvarMatrix(r$sol$itr$barx[[1]], prob$bardim[1])) );
print( list(X2=getbarvarMatrix(r$sol$itr$barx[[2]], prob$bardim[2])) );
}
sdo2();
sdo_lmi.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : sdo_lmi.R
#
# 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
##
library("Rmosek")
getbarvarMatrix <- function(barvar, bardim)
{
N <- as.integer(bardim)
new("dspMatrix", x=barvar, uplo="L", Dim=c(N,N))
}
sdo_lmi <- function()
{
# Specify the non-matrix variable part of the problem.
prob <- list(sense="min")
prob$c <- c(1, 1)
prob$c0 <- 1
# Specify variable bounds
prob$bx <- rbind(blx=rep(-Inf,2), bux=rep( Inf,2))
# The following two entries must always be defined, even if set to zero.
prob$A <- Matrix(c(0, 0), nrow=1, sparse=TRUE)
prob$bc <- rbind(blc=rep(-Inf,1),buc=rep(Inf, 1))
prob$F <- sparseMatrix(i=c(1, 1, 2, 3, 3, 4),
j=c(1, 2, 2, 1, 2, 1),
x=c(-1, -1, 3, sqrt(2), sqrt(2), 3), dims=c(4, 2))
prob$g <- c(0, -1, 0, -1)
prob$cones <- matrix(list(), nrow=3, ncol=2)
prob$cones[,1] <- list("RPLUS", 1, NULL)
prob$cones[,2] <- list("SVEC_PSD_CONE", 3, NULL)
# Specify semidefinite matrix variables (one 2x2 block)
prob$bardim <- c(2)
# Block triplet format specifying the lower triangular part
# of the symmetric coefficient matrix 'barc':
prob$barc$j <- c(1, 1, 1)
prob$barc$k <- c(1, 2, 2)
prob$barc$l <- c(1, 1, 2)
prob$barc$v <- c(1, 0, 1)
# Block triplet format specifying the lower triangular part
# of the symmetric coefficient matrix 'barF' for the ACC:
prob$barF$i <- c(1, 1)
prob$barF$j <- c(1, 1)
prob$barF$k <- c(1, 2)
prob$barF$l <- c(1, 1)
prob$barF$v <- c(0, 1)
# Solve the problem
r <- mosek(prob)
# Print matrix variable and return the solution
stopifnot(identical(r$response$code, 0))
print( list(barx=getbarvarMatrix(r$sol$itr$barx[[1]], prob$bardim[1]), xx=r$sol$itr$xx) )
r$sol
}
sdo_lmi()
simple.R
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : simple.R
#
# Purpose : To demonstrate a very simple example using the Rmosek package
# to read a problem file, solve the optimization problem and
# write the solution.
##
library("Rmosek")
simple <- function(filename)
{
# Read problem
r <- mosek_read(filename, list(verbose=1, usesol=FALSE, useparam=TRUE))
stopifnot(identical(r$response$code, 0))
prob <- r$prob
# Solve problem
r <- mosek(prob, list(verbose=1))
stopifnot(identical(r$response$code, 0))
sol <- r$sol
# Print solution
print(sol)
}
# Run only if called directly from the command-line
if( sys.nframe() == 0L )
{
args <- commandArgs(trailingOnly=TRUE)
if( !file.exists(args[1]) ) {
stop("Expected problem file as input argument.")
}
simple(args[1])
}
solutionquality.R
###
## Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
##
## File : solutionquality.R
##
## Purpose : To demonstrate how to examine the quality of a solution.
###
library("Rmosek")
solutionquality <- function()
{
# Specify problem.
prob <- list(sense="max")
prob$c <- c(3, 1, 5, 1)
prob$A <- sparseMatrix(i=c(1, 1, 1, 2, 2, 2, 2, 3, 3),
j=c(1, 2, 3, 1, 2, 3, 4, 2, 4),
x=c(3, 1, 2, 2, 1, 3, 1, 2, 3))
prob$bc <- rbind(blc=c(30, 15, -Inf),
buc=c(30, Inf, 25))
prob$bx <- rbind(blx=rep(0,4),
bux=c(Inf, 10, Inf, Inf))
# Solve the problem.
r <- mosek(prob, list(verbose=1, soldetail=2))
stopifnot(identical(r$response$code, 0))
# Examine interior-point solution quality.
isol <- r$sol$itr
if( is.element(isol$solsta, list("OPTIMAL")) )
{
pobj <- isol$pobjval
dobj <- isol$dobjval
abs_obj_gap <- abs(dobj-pobj)
rel_obj_gap <- abs_obj_gap/(1.0 + min(abs(pobj),abs(dobj)))
max_primal_viol <- max( isol$maxinfeas$pcon,
isol$maxinfeas$pbound,
isol$maxinfeas$pbarvar,
isol$maxinfeas$pcone )
max_dual_viol <- max( isol$maxinfeas$dcon,
isol$maxinfeas$dbound,
isol$maxinfeas$dbarvar,
isol$maxinfeas$dcone )
## Assume the application needs the solution to be within
## 1e-6 of optimality in an absolute sense. Another approach
## would be looking at the relative objective gap
cat("Customized solution information.\n")
cat(" Absolute objective gap: ",abs_obj_gap,"\n")
cat(" Relative objective gap: ",rel_obj_gap,"\n")
cat(" Max primal violation : ",max_primal_viol,"\n")
cat(" Max dual violation : ",max_dual_viol,"\n")
accepted <- TRUE
if( rel_obj_gap > 1e-6 ) {
print("Warning: The relative objective gap is LARGE.")
accepted <- FALSE
}
## We will accept a primal and dual infeasibility of 1e-6.
## These numbers should chosen problem dependent.
if( max_primal_viol > 1e-6 ) {
print("Warning: Primal violation is too LARGE")
accepted <- FALSE
}
if( max_dual_viol > 1e-6 ) {
print("Warning: Dual violation is too LARGE.")
accepted <- FALSE
}
if( accepted ) {
print("Optimal primal interior-point solution")
print(isol$xx)
}
} else if( is.element(isol$solsta, list("DUAL_INFEASIBLE_CER","PRIMAL_INFEASIBLE_CER")) ) {
cat("Primal or dual infeasibility certificate found.\n")
} else {
cat("Other solution status.\n")
}
}
solutionquality()
