17 List of examples¶
List of examples shipped in the distribution of Optimizer API for Python:
File |
Description |
|---|---|
A simple problem with one affine conic constraint (ACC) |
|
A simple problem with two affine conic constraints (ACC) |
|
Demonstrates the MOSEK interface to BLAS/LAPACK linear algebra routines |
|
An example of data/progress callback |
|
A simple conic exponential problem |
|
Implementation of a concurrent optimizer for linear and mixed-integer problems |
|
A simple conic quadratic problem |
|
A simple problem with disjunctive constraints (DJC) |
|
A simple example of how to repair an infeasible problem |
|
A simple geometric program (GP) in conic form |
|
A Hello World example |
|
A simple linear problem |
|
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 |
|
Uses MOSEK OptServer to solve an optimization problem asynchronously |
|
Uses MOSEK OptServer to solve an optimization problem synchronously |
|
Demonstrates parallel optimization using a batch method in MOSEK |
|
Shows how to set optimizer parameters and read information items |
|
Shows how to obtain and analyze a primal infeasibility certificate |
|
Portfolio optimization - basic Markowitz model |
|
Portfolio optimization - efficient frontier |
|
Portfolio optimization - market impact costs |
|
Portfolio optimization - transaction costs |
|
Portfolio optimization - cardinality constraints |
|
Portfolio optimization - factor model |
|
A simple power cone problem |
|
A simple quadratically constrained quadratic problem |
|
A simple quadratic problem |
|
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. |
|
Sensitivity analysis performed on a small linear problem |
|
A simple I/O example: read problem from a file, solve and write solutions |
|
Demonstrates how to examine the quality of a solution |
|
Demonstrates solving a linear system with the basis matrix |
|
Demonstrates solving a general linear system |
|
Shows how to find a Cholesky factorization of a sparse matrix |
Additional examples can be found on the MOSEK website and in other MOSEK publications.
acc1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : acc1.py
#
# Purpose : Tutorial example for affine conic constraints.
# Models the problem:
#
# maximize c^T x
# subject to sum(x) = 1
# gamma >= |Gx+h|_2
##
import sys, mosek
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
# Define problem data
n, k = 3, 2
# Only a symbolic constant
inf = 0.0
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
# Create n free variables
task.appendvars(n)
task.putvarboundsliceconst(0, n, mosek.boundkey.fr, -inf, inf)
# Set up the objective
c = [2, 3, -1]
task.putobjsense(mosek.objsense.maximize)
task.putclist(range(n), c)
# One linear constraint - sum(x) = 1
task.appendcons(1)
task.putarow(0, range(n), [1]*n)
task.putconbound(0, mosek.boundkey.fx, 1.0, 1.0)
# Append empty AFE rows for affine expression storage
task.appendafes(k + 1)
# G matrix in sparse form
Gsubi = [0, 0, 1, 1]
Gsubj = [0, 1, 0, 2]
Gval = [1.5, 0.1, 0.3, 2.1]
# Other data
h = [0, 0.1]
gamma = 0.03
# Construct F matrix in sparse form
Fsubi = [i + 1 for i in Gsubi] # G will be placed from row number 1 in F
Fsubj = Gsubj
Fval = Gval
# Fill in F storage
task.putafefentrylist(Fsubi, Fsubj, Fval)
# Fill in g storage
task.putafeg(0, gamma)
task.putafegslice(1, k+1, h)
# Define a conic quadratic domain
quadDom = task.appendquadraticconedomain(k + 1)
# Create the ACC
task.appendacc(quadDom, # Domain index
range(k+1), # Indices of AFE rows [0,...,k]
None) # Ignored
task.writedata("acc1.ptf")
# Solve and retrieve solution
task.optimize()
xx = task.getxx(mosek.soltype.itr)
if task.getsolsta(mosek.soltype.itr) == mosek.solsta.optimal:
print("Solution: {xx}".format(xx=list(xx)))
# Demonstrate retrieving activity of ACC
activity = task.evaluateacc(mosek.soltype.itr,
0) # ACC index
print("Activity of ACC:: {activity}".format(activity=list(activity)))
# Demonstrate retrieving the dual of ACC
doty = task.getaccdoty(mosek.soltype.itr,
0) # ACC index
print("Dual of ACC:: {doty}".format(doty=list(doty)))
acc2.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : acc2.py
#
# Purpose : Tutorial example for affine conic constraints.
# Models the problem:
#
# maximize c^T x
# subject to sum(x) = 1
# gamma >= |Gx+h|_2
#
# This version inputs the linear constraint as an affine conic constraint.
##
import sys, mosek
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
# Define problem data
n, k = 3, 2
# Only a symbolic constant
inf = 0.0
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
# Create n free variables
task.appendvars(n)
task.putvarboundsliceconst(0, n, mosek.boundkey.fr, -inf, inf)
# Set up the objective
c = [2, 3, -1]
task.putobjsense(mosek.objsense.maximize)
task.putclist(range(n), c)
# Set AFE rows representing the linear constraint
task.appendafes(1)
task.putafefrow(0, range(n), [1]*n)
task.putafeg(0, -1.0)
# Set AFE rows representing the quadratic constraint
task.appendafes(k + 1)
task.putafefrow(2, # afeidx, row number
[0, 1], # varidx, column numbers
[1.5, 0.1]) # values
task.putafefrow(3, # afeidx, row number
[0, 2], # varidx, column numbers
[0.3, 2.1]) # values
h = [0, 0.1]
gamma = 0.03
task.putafeg(1, gamma)
task.putafegslice(2, k+2, h)
# Define domains
zeroDom = task.appendrzerodomain(1)
quadDom = task.appendquadraticconedomain(k + 1)
# Append affine conic constraints
task.appendacc(zeroDom, # Domain index
[0], # Indices of AFE rows
None) # Ignored
task.appendacc(quadDom, # Domain index
[1,2,3], # Indices of AFE rows
None) # Ignored
# Solve and retrieve solution
task.optimize()
xx = task.getxx(mosek.soltype.itr)
if task.getsolsta(mosek.soltype.itr) == mosek.solsta.optimal:
print("Solution: {xx}".format(xx=list(xx)))
# Demonstrate retrieving activity of ACC
activity = task.evaluateacc(mosek.soltype.itr,
1) # ACC index
print("Activity of quadratic ACC:: {activity}".format(activity=list(activity)))
# Demonstrate retrieving the dual of ACC
doty = task.getaccdoty(mosek.soltype.itr,
1) # ACC index
print("Dual of quadratic ACC:: {doty}".format(doty=list(doty)))
blas_lapack.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : blas_lapack.py
#
# Purpose : To demonstrate how to call BLAS/LAPACK routines
# for whose MOSEK provides simplified interfaces.
##
import mosek
def print_matrix(x, r, c):
for i in range(r):
print([x[j * r + i] for j in range(c)])
with mosek.Env() as env:
n = 3
m = 2
k = 3
alpha = 2.0
beta = 0.5
x = [1.0, 1.0, 1.0]
y = [1.0, 2.0, 3.0]
z = [1.0, 1.0]
v = [0.0, 0.0]
#A has m=2 rows and k=3 cols
A = [1.0, 1.0, 2.0, 2.0, 3., 3.]
#B has k=3 rows and n=3 cols
B = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
C = [0.0 for i in range(n * m)]
D = [1.0, 1.0, 1.0, 1.0]
Q = [1.0, 0.0, 0.0, 2.0]
# BLAS routines
xy = env.dot(n, x, y)
print("dot results= %f\n" % xy)
env.axpy(n, alpha, x, y)
print("\naxpy results is ")
print_matrix(y, 1, len(y))
env.gemv(mosek.transpose.no, m, n, alpha, A, x, beta, z)
print("\ngemv results is ")
print_matrix(z, 1, len(z))
env.gemm(mosek.transpose.no, mosek.transpose.no,
m, n, k, alpha, A, B, beta, C)
print("\ngemm results is ")
print_matrix(C, m, n)
env.syrk(mosek.uplo.lo, mosek.transpose.no, m, k, alpha, A, beta, D)
print("\nsyrk results is")
print_matrix(D, m, m)
# LAPACK routines
env.potrf(mosek.uplo.lo, m, Q)
print("\npotrf results is ")
print_matrix(Q, m, m)
env.syeig(mosek.uplo.lo, m, Q, v)
print("\nsyeig results is")
print_matrix(v, 1, m)
env.syevd(mosek.uplo.lo, m, Q, v)
print("\nsyevd results is")
print('v: ')
print_matrix(v, 1, m)
print('Q: ')
print_matrix(Q, m, m)
print("Exiting...")
callback.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : callback.py
#
# Purpose : To demonstrate how to use the progress
# callback.
#
# Use this script as follows:
#
# callback.py psim 25fv47.mps
# callback.py dsim 25fv47.mps
# callback.py intpnt 25fv47.mps
#
# The first argument tells which optimizer to use
# i.e. psim is primal simplex, dsim is dual simplex
# and intpnt is interior-point.
##
import sys
import numpy
import mosek
from mosek import *
def makeUserCallback(maxtime, task):
pass
def userCallback(caller,
douinf,
intinf,
lintinf):
opttime = 0.0
if caller == callbackcode.begin_intpnt:
print("Starting interior-point optimizer")
elif caller == callbackcode.intpnt:
itrn = intinf[iinfitem.intpnt_iter]
pobj = douinf[dinfitem.intpnt_primal_obj]
dobj = douinf[dinfitem.intpnt_dual_obj]
stime = douinf[dinfitem.intpnt_time]
opttime = douinf[dinfitem.optimizer_time]
print("Iterations: %-3d" % itrn)
print(" Elapsed time: %6.2f(%.2f) " % (opttime, stime))
print(" Primal obj.: %-18.6e Dual obj.: %-18.6e" % (pobj, dobj))
elif caller == callbackcode.end_intpnt:
print("Interior-point optimizer finished.")
elif caller == callbackcode.begin_primal_simplex:
print("Primal simplex optimizer started.")
elif caller == callbackcode.update_primal_simplex:
itrn = intinf[iinfitem.sim_primal_iter]
pobj = douinf[dinfitem.sim_obj]
stime = douinf[dinfitem.sim_time]
opttime = douinf[dinfitem.optimizer_time]
print("Iterations: %-3d" % itrn)
print(" Elapsed time: %6.2f(%.2f)" % (opttime, stime))
print(" Obj.: %-18.6e" % pobj)
elif caller == callbackcode.end_primal_simplex:
print("Primal simplex optimizer finished.")
elif caller == callbackcode.begin_dual_simplex:
print("Dual simplex optimizer started.")
elif caller == callbackcode.update_dual_simplex:
itrn = intinf[iinfitem.sim_dual_iter]
pobj = douinf[dinfitem.sim_obj]
stime = douinf[dinfitem.sim_time]
opttime = douinf[dinfitem.optimizer_time]
print("Iterations: %-3d" % itrn)
print(" Elapsed time: %6.2f(%.2f)" % (opttime, stime))
print(" Obj.: %-18.6e" % pobj)
elif caller == callbackcode.end_dual_simplex:
print("Dual simplex optimizer finished.")
elif caller == callbackcode.new_int_mio:
print("New integer solution has been located.")
xx = task.getxx(soltype.itg)
print(xx)
print("Obj.: %f" % douinf[dinfitem.mio_obj_int])
else:
pass
if opttime >= maxtime:
# mosek is spending too much time. Terminate it.
print("Terminating.")
return 1
return 0
return userCallback
def main(args):
# To run a continuous problem example
filename = "../data/25fv47.mps"
slvr = "intpnt"
# To run a mixed-integer example
#filename = "../data/milo1.lp"
#slvr = ""
if len(args) < 3:
print("Usage: callback ( psim | dsim | intpnt ) filename")
if len(args) > 1:
slvr = args[1]
if len(args) > 2:
filename = args[2]
with mosek.Task() as task:
task.readdata(filename)
if slvr == 'psim':
task.putintparam(iparam.optimizer,
optimizertype.primal_simplex)
elif slvr == "dsim":
task.putintparam(iparam.optimizer, optimizertype.dual_simplex)
elif slvr == "intpnt":
task.putintparam(iparam.optimizer, optimizertype.intpnt)
usercallback = makeUserCallback(maxtime=0.05, task=task)
task.set_InfoCallback(usercallback)
task.optimize()
if __name__ == '__main__':
main(sys.argv)
ceo1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : ceo1.py
#
# Purpose : Demonstrates how to solve small conic exponential
# optimization problem using the MOSEK Python API.
##
import sys
import mosek
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Only a symbolic constant
inf = 0.0
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
c = [1.0, 1.0, 0.0]
a = [1.0, 1.0, 1.0]
numvar, numcon = 3, 1
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
# Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
# Set up the linear part of the problem
task.putcslice(0, numvar, c)
task.putarow(0, [0, 1, 2], a)
task.putvarboundslice(0, numvar, [mosek.boundkey.fr] * numvar, [inf] * numvar, [inf] * numvar)
task.putconbound(0, mosek.boundkey.fx, 1.0, 1.0)
# Add a conic constraint
# Create a 3x3 identity matrix F
task.appendafes(3)
task.putafefentrylist([0, 1, 2], # Rows
[0, 1, 2], # Columns
[1.0] * 3)
# Exponential cone (x(0),x(1),x(2)) \in EXP
expdomain = task.appendprimalexpconedomain()
task.appendacc(expdomain, # Domain
[0, 1, 2], # Rows from F
None) # Unused
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.minimize)
task.writedata("ceo1-py.ptf")
# Optimize the task
task.optimize()
# Print a summary containing information
# about the solution for debugging purposes
task.solutionsummary(mosek.streamtype.msg)
prosta = task.getprosta(mosek.soltype.itr)
solsta = task.getsolsta(mosek.soltype.itr)
# Output a solution
xx = task.getxx(mosek.soltype.itr)
if solsta == mosek.solsta.optimal:
print("Optimal solution: %s" % xx)
elif solsta == mosek.solsta.dual_infeas_cer:
print("Primal or dual infeasibility.\n")
elif solsta == mosek.solsta.prim_infeas_cer:
print("Primal or dual infeasibility.\n")
elif mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.code))
if msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
concurrent1.py
#
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: concurrent1.py
#
# Purpose: Demonstrates a simple implementation of a concurrent optimizer.
#
# The concurrent optimizer starts a few parallel optimizations
# of the same problem using different algorithms, and reports
# a solution when the first optimizer is ready.
#
# This example also demonstrates how to define a simple callback handler
# that stops the optimizer when requested.
#
import mosek, sys
from threading import Thread
# Defines a Mosek callback function whose only function
# is to indicate if the optimizer should be stopped.
stop = False
firstStop = -1
def cbFun(code):
return 1 if stop else 0
# Stream handling method
def streamwriter(s):
sys.stdout.write('{0}'.format(s))
# firstOK, res, trm = optimize(tasks)
#
# Takes a list of tasks and optimizes then in parallel. The
# response code and termination code from each optimization is
# returned in ``res`` and ``trm``.
#
# When one task completes with rescode == ok, others are terminated.
#
# firstOK is the index of the first task that returned
# with rescode == ok. Whether or not this task contains the
# most valuable answer, is for the caller to decide. If none
# completed without error returns -1.
def runTask(num, task, res, trm):
global stop
global firstStop
try:
trm[num] = task.optimize();
res[num] = mosek.rescode.ok
except mosek.MosekException as e:
trm[num] = mosek.rescode.err_unknown
res[num] = e.errno
finally:
# If this finished with success, inform other tasks to interrupt
if res[num] == mosek.rescode.ok:
if not stop:
firstStop = num
stop = True
def optimize(tasks):
n = len(tasks)
res = [ mosek.rescode.err_unknown ] * n
trm = [ mosek.rescode.err_unknown ] * n
# Set a callback function
for t in tasks:
t.set_Progress(cbFun)
# Start parallel optimizations, one per task
jobs = [ Thread(target=runTask, args=(i, tasks[i], res, trm)) for i in range(n) ]
for j in jobs:
j.start()
for j in jobs:
j.join()
# For debugging, print res and trm codes for all optimizers
for i in range(n):
print("Optimizer {0} res {1} trm {2}".format(i, res[i], trm[i]))
return firstStop, res, trm
#
# idx, winTask, winTrm, winRes = optimizeconcurrent(task, optimizers)
#
# Given a continuous task, set up jobs to optimize it
# with a list of different solvers.
#
# Returns an index, corresponding to the optimization
# task that is returned as winTask. This is the task
# with the best possible status of those that finished.
# If none task is considered successful returns -1.
def optimizeconcurrent(task, optimizers):
n = len(optimizers)
tasks = [ mosek.Task(task) for _ in range(n) ]
# Choose various optimizers for cloned tasks
for i in range(n):
tasks[i].putintparam(mosek.iparam.optimizer, optimizers[i])
# Solve tasks in parallel
firstOK, res, trm = optimize(tasks)
if firstOK >= 0:
return firstOK, tasks[firstOK], trm[firstOK], res[firstOK]
else:
return -1, None, None, None
#
# idx, winTask, winTrm, winRes = optimizeconcurrent(task, optimizers)
#
# Given a mixed-integer task, set up jobs to optimize it
# with different values of seed. That will lead to
# different execution paths of the optimizer.
#
# Returns an index, corresponding to the optimization
# task that is returned as winTask. This is the task
# with the best value of the objective function.
# If none task is considered successful returns -1.
#
# Typically, the input task would contain a time limit. The two
# major scenarios are:
# 1. Some clone ends before time limit - then it has optimum.
# 2. All clones reach time limit - pick the one with best objective.
def optimizeconcurrentMIO(task, seeds):
n = len(seeds)
tasks = [ mosek.Task(task) for _ in range(n) ]
# Choose various seeds for cloned tasks
for i in range(n):
tasks[i].putintparam(mosek.iparam.mio_seed, seeds[i])
# Solve tasks in parallel
firstOK, res, trm = optimize(tasks)
if firstOK >= 0:
# Pick the task that ended with res = ok
# and contains an integer solution with best objective value
sense = task.getobjsense();
bestObj = 1.0e+10 if sense == mosek.objsense.minimize else -1.0e+10
bestPos = -1
for i in range(n):
print("{0} {1}".format(i,tasks[i].getprimalobj(mosek.soltype.itg)))
for i in range(n):
if ((res[i] == mosek.rescode.ok) and
(tasks[i].getsolsta(mosek.soltype.itg) == mosek.solsta.prim_feas or
tasks[i].getsolsta(mosek.soltype.itg) == mosek.solsta.integer_optimal) and
((tasks[i].getprimalobj(mosek.soltype.itg) < bestObj)
if (sense == mosek.objsense.minimize) else
(tasks[i].getprimalobj(mosek.soltype.itg) > bestObj))):
bestObj = tasks[i].getprimalobj(mosek.soltype.itg)
bestPos = i
if bestPos >= 0:
return bestPos, tasks[bestPos], trm[bestPos], res[bestPos]
return -1, None, None, None
# This is an example of how one can use the methods
# optimizeconcurrent
# optimizeconcurrentMIO
#
# argv[0] : name of file with input problem
# argv[1]: (optional) time limit
def main(argv):
with mosek.Env() as env:
with mosek.Task(env, 0, 0) as task:
if len(argv) >= 2:
task.readdata(argv[1])
else:
task.readdata("../data/25fv47.mps")
# Optional time limit
if len(argv) >= 3:
task.putdouparam(mosek.dparam.optimizer_max_time, float(argv[2]))
if (task.getnumintvar() == 0):
# If the problem is continuous
# optimize it with three continuous optimizers.
# (Simplex will fail for non-linear problems)
optimizers = [
mosek.optimizertype.conic,
mosek.optimizertype.dual_simplex,
mosek.optimizertype.primal_simplex
]
idx, t, trm, res = optimizeconcurrent(task, optimizers)
else:
# Mixed-integer problem.
# Try various seeds.
seeds = [ 42, 13, 71749373 ]
idx, t, trm, res = optimizeconcurrentMIO(task, seeds)
# Check results and print the best answer
if idx >= 0:
print("Result from optimizer with index {0}: res {1} trm {2}".format(idx, res, trm))
t.set_Stream(mosek.streamtype.log, streamwriter)
t.optimizersummary(mosek.streamtype.log)
t.solutionsummary(mosek.streamtype.log);
else:
print("All optimizers failed.")
if __name__ == "__main__":
main(sys.argv)
cqo1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : cqo1.py
#
# Purpose : Demonstrates how to solve small conic quadratic
# optimization problem using the MOSEK Python API.
##
import sys
import mosek
# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
bkc = [mosek.boundkey.fx]
blc = [1.0]
buc = [1.0]
c = [0.0, 0.0, 0.0,
1.0, 1.0, 1.0]
bkx = [mosek.boundkey.lo, mosek.boundkey.lo, mosek.boundkey.lo,
mosek.boundkey.fr, mosek.boundkey.fr, mosek.boundkey.fr]
blx = [0.0, 0.0, 0.0,
-inf, -inf, -inf]
bux = [inf, inf, inf,
inf, inf, inf]
asub = [[0], [0], [0]]
aval = [[1.0], [1.0], [2.0]]
numvar = len(bkx)
numcon = len(bkc)
NUMANZ = 4
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
#Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
for j in range(numvar):
# Set the linear term c_j in the objective.
task.putcj(j, c[j])
# Set the bounds on variable j
# blx[j] <= x_j <= bux[j]
task.putvarbound(j, bkx[j], blx[j], bux[j])
for j in range(len(aval)):
# Input column j of A
task.putacol(j, # Variable (column) index.
# Row index of non-zeros in column j.
asub[j],
aval[j]) # Non-zero Values of column j.
for i in range(numcon):
task.putconbound(i, bkc[i], blc[i], buc[i])
# Input the affine conic constraints
# Create a matrix F such that F * x = [x(3),x(0),x(1),x(4),x(5),x(2)]
task.appendafes(6)
task.putafefentrylist(range(6), # Rows
[3, 0, 1, 4, 5, 2], # Columns
[1.0] * 6)
# Quadratic cone (x(3),x(0),x(1)) \in QUAD_3
quadcone = task.appendquadraticconedomain(3)
task.appendacc(quadcone, # Domain
[0, 1, 2], # Rows from F
None) # Unused
# Rotated quadratic cone (x(4),x(5),x(2)) \in RQUAD_3
rquadcone = task.appendrquadraticconedomain(3);
task.appendacc(rquadcone, # Domain
[3, 4, 5], # Rows from F
None) # Unused
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.minimize)
# Optimize the task
task.optimize()
task.writedata("cqo1.ptf")
# Print a summary containing information
# about the solution for debugging purposes
task.solutionsummary(mosek.streamtype.msg)
prosta = task.getprosta(mosek.soltype.itr)
solsta = task.getsolsta(mosek.soltype.itr)
# Output a solution
xx = task.getxx(mosek.soltype.itr)
if solsta == mosek.solsta.optimal:
print("Optimal solution: %s" % xx)
elif solsta == mosek.solsta.dual_infeas_cer:
print("Primal or dual infeasibility.\n")
elif solsta == mosek.solsta.prim_infeas_cer:
print("Primal or dual infeasibility.\n")
elif mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
djc1.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: djc1.py
#
# 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
##
import sys
from mosek import *
# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0
def main():
# Create a task object
with Task() as task:
# Append free variables
numvar = 4
task.appendvars(numvar)
task.putvarboundsliceconst(0, numvar, boundkey.fr, -inf, inf)
# The linear part: the linear constraint
task.appendcons(1)
task.putarow(0, range(numvar), [1] * numvar)
task.putconbound(0, boundkey.lo, -10.0, -10.0)
# The linear part: objective
task.putobjsense(objsense.minimize)
task.putclist(range(numvar), [2, 1, 3, 1])
# Fill in the affine expression storage F, g
numafe = 10
task.appendafes(numafe)
fafeidx = [0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9]
fvaridx = [0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]
fval = [1.0, -2.0, 1.0, -3.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
g = [1.0, 2.0, 0.0, 0.0, 0.0, 0.0, -2.5, -2.5, -2.5, -2.5]
task.putafefentrylist(fafeidx, fvaridx, fval)
task.putafegslice(0, numafe, g)
# Create domains
zero1 = task.appendrzerodomain(1)
zero2 = task.appendrzerodomain(2)
rminus1 = task.appendrminusdomain(1)
# Append disjunctive constraints
numdjc = 2
task.appenddjcs(numdjc)
# First disjunctive constraint
task.putdjc(0, # DJC index
[rminus1, zero2, rminus1, zero2], # Domains (domidxlist)
[0, 4, 5, 1, 2, 3], # AFE indices (afeidxlist)
None, # Unused
[2, 2] ) # Term sizes (termsizelist)
# Second disjunctive constraint
task.putdjc(1, # DJC index
[zero1, zero1, zero1, zero1], # Domains (domidxlist)
[6, 7, 8, 9], # AFE indices (afeidxlist)
None, # Unused
[1, 1, 1, 1] ) # Term sizes (termidxlist)
# Useful for debugging
task.writedata("djc.ptf") # Write file in human-readable format
task.set_Stream(streamtype.log, sys.stdout.write) # Attach a log stream printer to the task
# Solve the problem
task.optimize()
# Print a summary containing information
# about the solution for debugging purposes
task.solutionsummary(streamtype.msg)
# Get status information about the solution
sta = task.getsolsta(soltype.itg)
if (sta == solsta.integer_optimal):
xx = task.getxx(soltype.itg)
print("Optimal solution: ")
for i in range(numvar):
print("x[" + str(i) + "]=" + str(xx[i]))
else:
print("Another solution status")
# call the main function
try:
main()
except Error as e:
print("ERROR: %s" % str(e.errno))
if e.msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
feasrepairex1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : feasrepairex1.py
#
# Purpose : To demonstrate how to use the MSK_relaxprimal function to
# locate the cause of an infeasibility.
#
# Syntax : On command line
# python feasrepairex1.py feasrepair.lp
# feasrepair.lp is located in mosek\<version>\tools\examples.
##
import sys
import mosek
# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main(inputfile):
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
# Read data
task.readdata(inputfile)
task.putintparam(mosek.iparam.log_feas_repair, 3)
task.primalrepair(None, None, None, None)
sum_viol = task.getdouinf(mosek.dinfitem.primal_repair_penalty_obj)
print("Minimized sum of violations = %e" % sum_viol)
task.optimize()
task.solutionsummary(mosek.streamtype.msg)
# call the main function
try:
filename = "../data/feasrepair.lp"
if len(sys.argv) > 1:
filename = sys.argv[1]
main(filename)
except Exception as e:
print(e)
raise
gp1.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: gp1.py
#
# 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
#
from numpy import log, exp
from mosek import *
import sys
# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
# 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 )
def max_volume_box(Aw, Af, alpha, beta, gamma, delta):
# Basic dimensions of our problem
numvar = 3 # Variables in original problem
x, y, z = 0, 1, 2 # Indices of variables
numcon = 3 # Linear constraints in original problem
# Linear part of the problem
cval = [1, 1, 1]
asubi = [0, 0, 1, 1, 2, 2]
asubj = [y, z, x, y, z, y]
aval = [1.0, 1.0, 1.0, -1.0, 1.0, -1.0]
bkc = [boundkey.up, boundkey.ra, boundkey.ra]
blc = [-inf, log(alpha), log(gamma)]
buc = [log(Af), log(beta), log(delta)]
with Task() as task:
task.set_Stream(streamtype.log, streamprinter)
# Add variables and constraints
task.appendvars(numvar)
task.appendcons(numcon)
# Objective is the sum of three first variables
task.putobjsense(objsense.maximize)
task.putcslice(0, numvar, cval)
task.putvarboundsliceconst(0, numvar, boundkey.fr, -inf, inf)
# Add the three linear constraints
task.putaijlist(asubi, asubj, aval)
task.putconboundslice(0, numvar, bkc, blc, buc)
# Affine expressions appearing in affine conic constraints
# in this order:
# u1, u2, x+y+log(2/Awall), x+z+log(2/Awall), 1.0, u1+u2-1.0
numafe = 6
u1, u2 = 3, 4 # Indices of slack variables
afeidx = [0, 1, 2, 2, 3, 3, 5, 5]
varidx = [u1, u2, x, y, x, z, u1, u2]
fval = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
gfull = [0, 0, log(2/Aw), log(2/Aw), 1.0, -1.0]
# New variables u1, u2
task.appendvars(2)
task.putvarboundsliceconst(u1, u2+1, boundkey.fr, -inf, inf)
# Append affine expressions
task.appendafes(numafe)
task.putafefentrylist(afeidx, varidx, fval)
task.putafegslice(0, numafe, gfull)
# Two affine conic constraints
expdom = task.appendprimalexpconedomain()
# (u1, 1, x+y+log(2/Awall)) \in EXP
task.appendacc(expdom, [0, 4, 2], None)
# (u2, 1, x+z+log(2/Awall)) \in EXP
task.appendacc(expdom, [1, 4, 3], None)
# The constraint u1+u2-1 \in \ZERO is added also as an ACC
task.appendacc(task.appendrzerodomain(1), [5], None)
# Solve and map to original h, w, d
task.optimize()
task.writedata("gp1.ptf");
xyz = task.getxxslice(soltype.itr, 0, numvar)
return exp(xyz)
Aw, Af, alpha, beta, gamma, delta = 200.0, 50.0, 2.0, 10.0, 2.0, 10.0
h,w,d = max_volume_box(Aw, Af, alpha, beta, gamma, delta)
print("h={0:.3f}, w={1:.3f}, d={2:.3f}".format(h, w, d))
helloworld.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: helloworld.py
#
# The most basic example of how to get started with MOSEK.
from mosek import *
with Task() as task: # Create Task
task.appendvars(1) # 1 variable x
task.putcj(0, 1.0) # c_0 = 1.0
task.putvarbound(0, boundkey.ra, 2.0, 3.0) # 2.0 <= x <= 3.0
task.putobjsense(objsense.minimize) # minimize
task.optimize() # Optimize
x = task.getxx(soltype.itr) # Get solution
print("Solution x = {}".format(x[0])) # Print solution
lo1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : lo1.py
#
# Purpose : Demonstrates how to solve small linear
# optimization problem using the MOSEK Python API.
##
import sys
import mosek
# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task object
with mosek.Task() as task:
# Attach a log stream printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
# Bound keys for constraints
bkc = [mosek.boundkey.fx,
mosek.boundkey.lo,
mosek.boundkey.up]
# Bound values for constraints
blc = [30.0, 15.0, -inf]
buc = [30.0, +inf, 25.0]
# Bound keys for variables
bkx = [mosek.boundkey.lo,
mosek.boundkey.ra,
mosek.boundkey.lo,
mosek.boundkey.lo]
# Bound values for variables
blx = [0.0, 0.0, 0.0, 0.0]
bux = [+inf, 10.0, +inf, +inf]
# Objective coefficients
c = [3.0, 1.0, 5.0, 1.0]
# Below is the sparse representation of the A
# matrix stored by column.
asub = [[0, 1],
[0, 1, 2],
[0, 1],
[1, 2]]
aval = [[3.0, 2.0],
[1.0, 1.0, 2.0],
[2.0, 3.0],
[1.0, 3.0]]
numvar = len(bkx)
numcon = len(bkc)
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
# Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
for j in range(numvar):
# Set the linear term c_j in the objective.
task.putcj(j, c[j])
# Set the bounds on variable j
# blx[j] <= x_j <= bux[j]
task.putvarbound(j, bkx[j], blx[j], bux[j])
# Input column j of A
task.putacol(j, # Variable (column) index.
asub[j], # Row index of non-zeros in column j.
aval[j]) # Non-zero Values of column j.
# Set the bounds on constraints.
# blc[i] <= constraint_i <= buc[i]
for i in range(numcon):
task.putconbound(i, bkc[i], blc[i], buc[i])
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.maximize)
# Solve the problem
task.optimize()
# Print a summary containing information
# about the solution for debugging purposes
task.solutionsummary(mosek.streamtype.msg)
# Get status information about the solution
solsta = task.getsolsta(mosek.soltype.bas)
if (solsta == mosek.solsta.optimal):
xx = task.getxx(mosek.soltype.bas)
print("Optimal solution: ")
for i in range(numvar):
print("x[" + str(i) + "]=" + str(xx[i]))
elif (solsta == mosek.solsta.dual_infeas_cer or
solsta == mosek.solsta.prim_infeas_cer):
print("Primal or dual infeasibility certificate found.\n")
elif solsta == mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.Error as e:
print("ERROR: %s" % str(e.errno))
if e.msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
lo2.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : lo2.py
#
# Purpose : Demonstrates how to solve small linear
# optimization problem using the MOSEK Python API.
##
import sys
import mosek
# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
# Bound keys for constraints
bkc = [mosek.boundkey.fx,
mosek.boundkey.lo,
mosek.boundkey.up]
# Bound values for constraints
blc = [30.0, 15.0, -inf]
buc = [30.0, +inf, 25.0]
# Bound keys for variables
bkx = [mosek.boundkey.lo,
mosek.boundkey.ra,
mosek.boundkey.lo,
mosek.boundkey.lo]
# Bound values for variables
blx = [0.0, 0.0, 0.0, 0.0]
bux = [+inf, 10.0, +inf, +inf]
# Objective coefficients
c = [3.0, 1.0, 5.0, 1.0]
# We input the A matrix column-wise
# asub contains row indexes
asub = [[0, 1, 2],
[0, 1, 2, 3],
[0, 3]]
# acof contains coefficients
aval = [[3.0, 1.0, 2.0],
[2.0, 1.0, 3.0, 1.0],
[2.0, 3.0]]
numvar = len(bkx)
numcon = len(bkc)
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
#Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
for j in range(numvar):
# Set the linear term c_j in the objective.
task.putcj(j, c[j])
# Set the bounds on variable j
# blx[j] <= x_j <= bux[j]
task.putvarbound(j, bkx[j], blx[j], bux[j])
for i in range(numcon):
task.putconbound(i, bkc[i], blc[i], buc[i])
# Input row i of A
task.putarow(i, # Row index.
# Column indexes of non-zeros in row i.
asub[i],
aval[i]) # Non-zero Values of row i.
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.maximize)
# Optimize the task
task.optimize()
# Print a summary containing information
#about the solution for debugging purposes
task.solutionsummary(mosek.streamtype.msg)
prosta = task.getprosta(mosek.soltype.bas)
solsta = task.getsolsta(mosek.soltype.bas)
# Output a solution
xx = task.getxx(mosek.soltype.bas)
if solsta == mosek.solsta.optimal:
print("Optimal solution: %s" % xx)
elif solsta == mosek.solsta.dual_infeas_cer:
print("Primal or dual infeasibility.\n")
elif solsta == mosek.solsta.prim_infeas_cer:
print("Primal or dual infeasibility.\n")
elif mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
if e.msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
logistic.py
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: logistic.py
#
# Purpose: Implements logistic regression with regulatization.
#
# Demonstrates using the exponential cone and log-sum-exp in Optimizer API.
#
# Plots an example for 2D datasets.
from mosek import *
import numpy as np
import sys, itertools
inf = 0.0
# Adds ACCs for t_i >= log ( 1 + exp((1-2*y[i]) * theta' * X[i]) )
# Adds auxiliary variables, AFE rows and constraints
def softplus(task, d, n, theta, t, X, y):
nvar = task.getnumvar()
ncon = task.getnumcon()
nafe = task.getnumafe()
task.appendvars(2*n) # z1, z2
task.appendcons(n) # z1 + z2 = 1
task.appendafes(4*n) #theta * X[i] - t[i], -t[i], z1[i], z2[i]
z1, z2 = nvar, nvar+n
zcon = ncon
thetaafe, tafe, z1afe, z2afe = nafe, nafe+n, nafe+2*n, nafe+3*n
for i in range(n):
task.putvarname(z1+i,f"z1[{i}]")
task.putvarname(z2+i,f"z2[{i}]")
# z1 + z2 = 1
task.putaijlist(range(zcon, zcon+n), range(z1, z1+n), [1]*n)
task.putaijlist(range(zcon, zcon+n), range(z2, z2+n), [1]*n)
task.putconboundsliceconst(zcon, zcon+n, boundkey.fx, 1, 1)
task.putvarboundsliceconst(nvar, nvar+2*n, boundkey.fr, -inf, inf)
# Affine conic expressions
afeidx, varidx, fval = [], [], []
## Thetas
for i in range(n):
for j in range(d):
afeidx.append(thetaafe + i)
varidx.append(theta + j)
fval.append(-X[i][j] if y[i]==1 else X[i][j])
# -t[i]
afeidx.extend([thetaafe + i for i in range(n)] + [tafe + i for i in range(n)])
varidx.extend([t + i for i in range(n)] + [t + i for i in range(n)])
fval.extend([-1.0]*(2*n))
# z1, z2
afeidx.extend([z1afe + i for i in range(n)] + [z2afe + i for i in range(n)])
varidx.extend([z1 + i for i in range(n)] + [z2 + i for i in range(n)])
fval.extend([1.0]*(2*n))
# Add the expressions
task.putafefentrylist(afeidx, varidx, fval)
# Add a single row with the constant expression "1.0"
oneafe = task.getnumafe()
task.appendafes(1)
task.putafeg(oneafe, 1.0)
# Add an exponential cone domain
expdomain = task.appendprimalexpconedomain()
# Conic constraints
acci = task.getnumacc()
for i in range(n):
task.appendacc(expdomain, [z1afe+i, oneafe, thetaafe+i], None)
task.appendacc(expdomain, [z2afe+i, oneafe, tafe+i], None)
task.putaccname(acci, f"z1:theta[{i}]")
task.putaccname(acci+1,f"z2:t[{i}]")
acci += 2
# Model logistic regression (regularized with full 2-norm of theta)
# X - n x d matrix of data points
# y - length n vector classifying training points
# lamb - regularization parameter
def logisticRegression(X, y, lamb=1.0):
n, d = int(X.shape[0]), int(X.shape[1]) # num samples, dimension
with Task() as task:
# Variables [r; theta; t; u]
nvar = 1+d+2*n
task.appendvars(nvar)
task.putvarboundsliceconst(0, nvar, boundkey.fr, -inf, inf)
r, theta, t = 0, 1, 1+d
task.putvarname(r,"r");
for j in range(d): task.putvarname(theta+j,f"theta[{j}]");
for j in range(n): task.putvarname(t+j,f"t[{j}]");
# Objective lambda*r + sum(t)
task.putobjsense(objsense.minimize)
task.putcj(r, lamb)
task.putclist(range(t, t+n), [1.0]*n)
# Softplus function constraints
softplus(task, d, n, theta, t, X, y);
# Regularization
# Append a sequence of linear expressions (r, theta) to F
numafe = task.getnumafe()
task.appendafes(1+d)
task.putafefentry(numafe, r, 1.0)
for i in range(d):
task.putafefentry(numafe + i + 1, theta + i, 1.0)
# Add the constraint
task.appendaccseq(task.appendquadraticconedomain(1+d), numafe, None)
# Solution
task.writedata('logistic.ptf')
task.optimize()
xx = task.getxxslice(soltype.itr, theta, theta+d)
return xx
# Map the 2d data through all monomials of degree <= d
def mapFeature(p,d):
return np.array([(p[0]**a)*(p[1]**b) for a,b in itertools.product(range(d+1), range(d+1)) if a+b<=d])
def mapFeatures(x,d):
return np.array([mapFeature(p,d) for p in x])
# Load the file and map using degree d monomials
# The file format is
# x_1, x_2, y
# for all training examples
def loaddata(filename):
# Read coordinates and y values
x, y = [], []
with open(filename, "r") as f:
for l in f.readlines():
num = l.split(',')
x.append([float(num[0]), float(num[1])])
y.append(int(num[2]))
return np.array(x), np.array(y)
# Plot some 2d results
def plot2d(x, y, d, theta):
import matplotlib
import matplotlib.pyplot as plt
pos = np.where(y==1)
neg = np.where(y==0)
plt.scatter(x[pos,0], x[pos,1], marker='o', color='b')
plt.scatter(x[neg,0], x[neg,1], marker='x', color='r')
u = np.linspace(-1, 1, 50)
v = np.linspace(-1, 1, 50)
z = np.zeros(shape=(len(u), len(v)))
for i in range(len(u)):
for j in range(len(v)):
z[i,j] = np.dot(mapFeature([u[i],v[j]], d), theta)
plt.contour(u, v, z.T, [0])
plt.show()
###############################################################################
# Example from documentation is contained in
# https://datascienceplus.com/wp-content/uploads/2017/02/ex2data2.txt
'''
for lamb in [1e-0,1e-2,1e-4]:
x, y = loaddata("ex2data2.txt")
X = mapFeatures(x, d=6)
theta = logisticRegression(env, X, y, lamb)
plot2d(x, y, 6, theta)
'''
# Example 2: discover a circle
x, y = [], []
for x1,x2 in itertools.product(np.linspace(-1, 1, 30),np.linspace(-1, 1, 30)):
x.append([x1,x2])
y.append(0 if x1**2+x2**2<0.69 else 1)
x, y = np.array(x), np.array(y)
X = mapFeatures(x, d=2)
# print(list(x[:,0]))
# print(list(x[:,1]))
# for i in range(len(x)):
# print(f"X[{i}] = {list(X[i,:])}")
theta = logisticRegression(X, y, 0.1)
print(theta)
mico1.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: mico1.py
#
# 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
##
import mosek, sys
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, streamprinter)
task.appendvars(3) # x, y, t
x, y, t = 0, 1, 2
task.putvarboundsliceconst(0, 3, mosek.boundkey.fr, -0.0, 0.0)
# Integrality constraints
task.putvartypelist([x,y], [mosek.variabletype.type_int]*2)
# Set up the affine expression
# x, x-3.8, y, t, 1.0
task.appendafes(5)
task.putafefentrylist([0, 1, 2, 3],
[x,x,y,t],
[1,1,1,1])
task.putafegslice(0, 5, [0, -3.8, 0, 0, 1.0])
# Add constraint (x-3.8, 1, y) \in \EXP
task.appendacc(task.appendprimalexpconedomain(), [1, 4, 2], None)
# Add constraint (t, x, y) \in \QUAD
task.appendacc(task.appendquadraticconedomain(3), [3, 0, 2], None)
# Objective
task.putobjsense(mosek.objsense.minimize)
task.putcj(t, 1)
# Optimize the task
task.optimize()
task.solutionsummary(mosek.streamtype.msg)
xx = task.getxxslice(mosek.soltype.itg, 0, 2)
print(xx)
milo1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : milo1.py
#
# Purpose : Demonstrates how to solve a small mixed
# integer linear optimization problem using the MOSEK Python API.
##
import sys
import mosek
# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
bkc = [mosek.boundkey.up, mosek.boundkey.lo]
blc = [-inf, -4.0]
buc = [250.0, inf]
bkx = [mosek.boundkey.lo, mosek.boundkey.lo]
blx = [0.0, 0.0]
bux = [inf, inf]
c = [1.0, 0.64]
asub = [[0, 1], [0, 1]]
aval = [[50.0, 3.0], [31.0, -2.0]]
numvar = len(bkx)
numcon = len(bkc)
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
#Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
for j in range(numvar):
# Set the linear term c_j in the objective.
task.putcj(j, c[j])
# Set the bounds on variable j
# blx[j] <= x_j <= bux[j]
task.putvarbound(j, bkx[j], blx[j], bux[j])
# Input column j of A
task.putacol(j, # Variable (column) index.
# Row index of non-zeros in column j.
asub[j],
aval[j]) # Non-zero Values of column j.
task.putconboundlist(range(numcon), bkc, blc, buc)
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.maximize)
# Define variables to be integers
task.putvartypelist([0, 1],
[mosek.variabletype.type_int,
mosek.variabletype.type_int])
# Set max solution time
task.putdouparam(mosek.dparam.mio_max_time, 60.0);
# Optimize the task
task.optimize()
task.writedata("milo1.ptf")
# Print a summary containing information
# about the solution for debugging purposes
task.solutionsummary(mosek.streamtype.msg)
prosta = task.getprosta(mosek.soltype.itg)
solsta = task.getsolsta(mosek.soltype.itg)
# Output a solution
xx = task.getxx(mosek.soltype.itg)
if solsta in [mosek.solsta.integer_optimal]:
print("Optimal solution: %s" % xx)
elif solsta == mosek.solsta.prim_feas:
print("Feasible solution: %s" % xx)
elif mosek.solsta.unknown:
if prosta == mosek.prosta.prim_infeas_or_unbounded:
print("Problem status Infeasible or unbounded.\n")
elif prosta == mosek.prosta.prim_infeas:
print("Problem status Infeasible.\n")
elif prosta == mosek.prosta.unkown:
print("Problem status unkown.\n")
else:
print("Other problem status.\n")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as msg:
#print "ERROR: %s" % str(code)
if msg is not None:
print("\t%s" % msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
mioinitsol.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : mioinitsol.py
#
# Purpose : Demonstrates how to solve a small mixed
# integer linear optimization problem using the MOSEK Python API.
##
import sys
import mosek
# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
bkc = [mosek.boundkey.up]
blc = [-inf, ]
buc = [2.5]
bkx = [mosek.boundkey.lo,
mosek.boundkey.lo,
mosek.boundkey.lo,
mosek.boundkey.lo]
blx = [0.0, 0.0, 0.0, 0.0]
bux = [inf, inf, inf, inf]
c = [7.0, 10.0, 1.0, 5.0]
asub = [0, 0, 0, 0]
acof = [1.0, 1.0, 1.0, 1.0]
ptrb = [0, 1, 2, 3]
ptre = [1, 2, 3, 4]
numvar = len(bkx)
numcon = len(bkc)
# Input linear data
task.inputdata(numcon, numvar,
c, 0.0,
ptrb, ptre, asub, acof,
bkc, blc, buc,
bkx, blx, bux)
# Input objective sense
task.putobjsense(mosek.objsense.maximize)
# Define variables to be integers
task.putvartypelist([0, 1, 2],
[mosek.variabletype.type_int]*3)
# Assign values to integer variables.
# (We only set a slice of xx)
task.putxxslice(mosek.soltype.itg, 0, 3, [1.0, 1.0, 0.0])
# Request constructing the solution from integer variable values
task.putintparam(mosek.iparam.mio_construct_sol, mosek.onoffkey.on)
# Optimize
task.optimize()
task.writedata("mioinitsol.ptf")
task.solutionsummary(mosek.streamtype.log)
if task.solutiondef(mosek.soltype.itg):
# Output a solution
xx = task.getxx(mosek.soltype.itg)
print("Integer optimal solution")
for j in range(0, numvar):
print("\tx[%d] = %e" % (j, xx[j]))
# Was the initial guess used?
constr = task.getintinf(mosek.iinfitem.mio_construct_solution)
constrVal = task.getdouinf(mosek.dinfitem.mio_construct_solution_obj)
print("Construct solution utilization: {0}\nConstruct solution objective: {1:.3f}\n".format(constr, constrVal))
else:
print("No integer solution is available.")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
if e.msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
opt_server_async.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : opt_server_async.py
#
# Purpose : Demonstrates how to use MOSEK OptServer
# to solve optimization problem asynchronously
##
import mosek
import sys
import time
def streamprinter(msg):
sys.stdout.write(msg)
sys.stdout.flush()
if len(sys.argv) < 4:
print("Missing argument, syntax is:")
print(" opt-server-async inputfile host:port numpolls [cert]")
else:
filename = sys.argv[1]
addr = sys.argv[2]
numpolls = int(sys.argv[3])
token = None
cert = None if len(sys.argv) < 5 else sys.argv[4]
with mosek.Task() as task:
print("reading task from file")
task.readdata(filename)
if cert is not None:
task.putstrparam(mosek.sparam.remote_tls_cert_path,cert)
print("Solve the problem remotely (async)")
token = task.asyncoptimize(addr,"")
print("Task token: %s" % token)
with mosek.Task() as task:
task.readdata(filename)
if cert is not None:
task.putstrparam(mosek.sparam.remote_tls_cert_path,cert)
task.set_Stream(mosek.streamtype.log, streamprinter)
i = 0
while i < numpolls:
time.sleep(0.1)
print("poll %d..." % i)
respavailable, res, trm = task.asyncpoll(addr,
"",
token)
print("done!")
if respavailable:
print("solution available!")
respavailable, res, trm = task.asyncgetresult(addr,
"",
token)
task.solutionsummary(mosek.streamtype.log)
break
i = i + 1
if i == numpolls:
print("max number of polls reached, stopping host.")
task.asyncstop(addr,"", token)
opt_server_sync.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : opt_server_sync.py
#
# Purpose : Demonstrates how to use MOSEK OptServer
# to solve optimization problem synchronously
##
import mosek
import sys
def streamprinter(msg):
sys.stdout.write(msg)
sys.stdout.flush()
if len(sys.argv) <= 2:
print("Missing argument, syntax is:")
print(" opt_server_sync inputfile addr [certpath]")
else:
inputfile = sys.argv[1]
serveraddr = sys.argv[2]
tlscert = None if len(sys.argv) < 4 else sys.argv[3]
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, streamprinter)
# We assume that a problem file was given as the first command
# line argument (received in `argv')
task.readdata(inputfile)
# Set OptServer URL
task.putoptserverhost(serveraddr)
# Path to certificate, if any
if tlscert is not None:
task.putstrparam(mosek.sparam.remote_tls_cert_path, tlscert)
# Solve the problem remotely, no access token
trm = task.optimize()
# Print a summary of the solution
task.solutionsummary(mosek.streamtype.log)
parallel.py
#
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: parallel.py
#
# Purpose: Demonstrates parallel optimization using optimizebatch()
#
import mosek, sys
# Example of how to use env.optimizebatch()
# Optimizes tasks whose names were read from command line.
def main(argv):
n = len(argv) - 1
tasks = []
threadpoolsize = 6 # Size of thread pool available for all tasks
with mosek.Env() as env:
# Set up some example list of tasks to optimize
for i in range(n):
t = mosek.Task(env, 0, 0)
t.readdata(argv[i+1])
# We can set the number of threads for each task
t.putintparam(mosek.iparam.num_threads, 2)
tasks.append(t)
# Optimize all the given tasks in parallel
trm, res = env.optimizebatch(False, # No race
-1.0, # No time limit
threadpoolsize,
tasks) # List of tasks to optimize
for i in range(n):
print("Task {0} res {1} trm {2} obj_val {3} time {4}".format(
i,
res[i],
trm[i],
tasks[i].getdouinf(mosek.dinfitem.intpnt_primal_obj),
tasks[i].getdouinf(mosek.dinfitem.optimizer_time)))
if __name__ == "__main__":
main(sys.argv)
parameters.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : parameters.py
#
# Purpose : Demonstrates a very simple example about how to get/set
# parameters with MOSEK Python API
##
import mosek
# Create the mosek task
with mosek.Task() as task:
print('Test MOSEK parameter get/set functions')
# Set log level (integer parameter)
task.putintparam(mosek.iparam.log, 1)
# Select interior-point optimizer... (integer parameter)
task.putintparam(mosek.iparam.optimizer, mosek.optimizertype.intpnt)
# ... without basis identification (integer parameter)
task.putintparam(mosek.iparam.intpnt_basis, mosek.basindtype.never)
# Set relative gap tolerance (double parameter)
task.putdouparam(mosek.dparam.intpnt_co_tol_rel_gap, 1.0e-7)
# The same using explicit string names
task.putparam ("MSK_DPAR_INTPNT_CO_TOL_REL_GAP", "1.0e-7")
task.putnadouparam("MSK_DPAR_INTPNT_CO_TOL_REL_GAP", 1.0e-7 )
# Incorrect value
try:
task.putdouparam(mosek.dparam.intpnt_co_tol_rel_gap, -1.0)
except:
print('Wrong parameter value')
param = task.getdouparam(mosek.dparam.intpnt_co_tol_rel_gap)
print('Current value for parameter intpnt_co_tol_rel_gap = {}'.format(param))
# Define and solve an optimization problem here
# task.optimize()
# After optimization:
print('Get MOSEK information items')
tm = task.getdouinf(mosek.dinfitem.optimizer_time)
it = task.getintinf(mosek.iinfitem.intpnt_iter)
print('Time: {0}\nIterations: {1}'.format(tm,it))
pinfeas.py
##
# File : pinfeas.py
#
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# Purpose: Demonstrates how to fetch a primal infeasibility certificate
# for a linear problem
##
from mosek import *
import sys
# Set up a simple linear problem from the manual for test purposes
def testProblem():
inf = 0.0
task = Task()
task.appendvars(7)
task.appendcons(7)
task.putclist(range(0, 7), [1, 2, 5, 2, 1, 2, 1])
task.putaijlist([0,0,1,1,2,2,2,3,3,4,5,5,6,6],
[0,1,2,3,4,5,6,0,4,1,2,5,3,6],
[1] * 14)
task.putconboundslice(0, 7, [boundkey.up]*3+[boundkey.fx]*4,
[-inf, -inf, -inf, 1100, 200, 500, 500],
[200, 1000, 1000, 1100, 200, 500, 500])
task.putvarboundsliceconst(0, 7, boundkey.lo, 0, +inf)
return task
# Analyzes and prints infeasibility contributing elements
# sl - dual values for lower bounds
# su - dual values for upper bounds
# eps - tolerance for when a nunzero dual value is significant
def analyzeCertificate(sl, su, eps):
n = len(sl)
for i in range(n):
if abs(sl[i]) > eps:
print(f"#{i}: lower, dual = {sl[i]}")
if abs(su[i]) > eps:
print(f"#{i}: upper, dual = {su[i]}")
def pinfeas():
# In this example we set up a simple problem
# One could use any task or a task read from a file
task = testProblem()
# Useful for debugging
task.writedata("pinfeas.ptf")
task.set_Stream(streamtype.log, sys.stdout.write)
# Perform the optimization.
task.optimize()
task.solutionsummary(streamtype.log)
# Check problem status, we use the interior point solution
if task.getprosta(soltype.itr) == prosta.prim_infeas:
# Set the tolerance at which we consider a dual value as essential
eps = 1e-7
print("Variable bounds important for infeasibility: ")
analyzeCertificate(task.getslx(soltype.itr), task.getsux(soltype.itr), eps)
print("Constraint bounds important for infeasibility: ")
analyzeCertificate(task.getslc(soltype.itr), task.getsuc(soltype.itr), eps)
else:
print("The problem is not primal infeasible, no certificate to show");
pinfeas()
portfolio_1_basic.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_1_basic.py
#
# Purpose : Implements a basic portfolio optimization model.
##
import mosek
import sys
import numpy as np
if __name__ == '__main__':
n = 8
w = 59.0
mu = [0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628]
x0 = [8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0]
gamma = 36.0
GT = [
[0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638],
[0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506],
[0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914],
[0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742],
[0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ],
[0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187],
[0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327],
[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ]
]
k = len(GT)
inf = 0.0 # This value has no significance
# Variable offsets
numvar = n
voff_x = 0
# Constraints offsets
numcon = 1
coff_bud = 0
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Holding variable x of length n
# No other auxiliary variables are needed in this formulation
task.appendvars(numvar)
# Optionally we can give the variables names
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (1 + j))
# No short-selling in this model, all of x >= 0
task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)
# One linear constraint: total budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
rtemp = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp) # equals w + sum(x0)
# Input (gamma, GTx) in the AFE (affine expression) storage
# We need k+1 rows
task.appendafes(k + 1)
# The first affine expression = gamma
task.putafeg(0, gamma)
# The remaining k expressions comprise GT*x, we add them row by row
# In more realisic scenarios it would be better to extract nonzeros and input in sparse form
for i in range(0, k):
task.putafefrow(i + 1, range(voff_x, voff_x + n), GT[i])
# Input the affine conic constraint (gamma, GT*x) \in QCone
# Add the quadratic domain of dimension k+1
qdom = task.appendquadraticconedomain(k + 1)
# Add the constraint
task.appendaccseq(qdom, 0, None)
task.putaccname(0, "risk")
# Objective: maximize expected return mu^T x
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
# Solve the problem
task.optimize()
# Display solution summary for quick inspection of results.
# In this simplified example we skip checks for problem and solution status
task.solutionsummary(mosek.streamtype.msg)
# Check if the interior point solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itr)
if (solsta != mosek.solsta.optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
# Retrieve results
xx = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
expret = task.getprimalobj(mosek.soltype.itr)
print(f'Expected return: {expret:.10e} Std. deviation: {gamma:.4e}')
np.set_printoptions(precision=4)
print(f'Optimal portfolio: {np.array(xx)}')
portfolio_2_frontier.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_2_frontier.py
#
# Purpose : Implements a basic portfolio optimization model.
# Computes points on the efficient frontier.
##
import mosek
import sys
import numpy as np
if __name__ == '__main__':
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 ]
]
k = len(GT)
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]
inf = 0.0 # This value has no significance
# Offset of variables into the API variable.
numvar = n + 1
voff_x = 0
voff_s = n
# Offset of constraints
coff_bud = 0
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Variables:
task.appendvars(numvar)
# Optionally we can give the variables names
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (1 + j))
task.putvarname(voff_s, "s")
# No short-selling in this model, all of x >= 0
task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)
# s is free variable
task.putvarbound(voff_s, mosek.boundkey.fr, -inf, inf)
# One linear constraint: total budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
rtemp = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp) # equals w + sum(x0)
# Input (gamma, GTx) in the AFE (affine expression) storage
# We build the following F and g for variables [x, s]:
# [0, 1] [0 ]
# F = [0, 0], g = [0.5]
# [GT,0] [0 ]
# We need k+2 rows
task.appendafes(k + 2)
# The first affine expression is variable s (last variable, index n)
task.putafefrow(0, [n], [1.0])
# The second affine expression is constant 0.5
task.putafeg(1, 0.5)
# The remaining k affine expressions comprise GT*x, we add them row by row
# In more realisic scenarios it would be better to extract nonzeros and input in sparse form
for i in range(0, k):
task.putafefrow(i + 2, range(voff_x, voff_x + n), GT[i])
# Input the affine conic constraint (s, 0.5, GT*x) \in RQCone
# Add the quadratic domain of dimension k+1
rqdom = task.appendrquadraticconedomain(k + 2)
# Add the constraint
task.appendaccseq(rqdom, 0, None)
task.putaccname(0, "risk")
# Set objective coefficients (x part): mu'x - alpha * s
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Turn all log output off.
task.putintparam(mosek.iparam.log, 0)
for alpha in alphas:
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
task.putcj(voff_s, -alpha)
task.optimize()
# Display the solution summary for quick inspection of results.
# task.solutionsummary(mosek.streamtype.msg)
# Check if the interior point solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itr)
if solsta in [mosek.solsta.optimal]:
expret = 0.0
x = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
for j in range(0, n):
expret += mu[j] * x[j]
stddev = np.sqrt(task.getxxslice(mosek.soltype.itr, voff_s, voff_s + 1))
print("alpha = {0:.2e} exp. ret. = {1:.3e}, std. dev. {2:.3e}".format(alpha, expret, stddev[0]))
else:
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
print("An error occurred when solving for alpha=%e\n" % alpha)
raise Exception(f"Unexpected solution status: {solsta}")
portfolio_3_impact.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_3_impact.py
#
# Purpose : Implements a basic portfolio optimization model
# with x^(3/2) market impact costs.
##
import mosek
import sys
if __name__ == '__main__':
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 ]
]
gamma = 0.36
k = len(GT)
m = n * [0.01]
# This value has no significance.
inf = 0.0
# Offset of variables.
numvar = 3 * n
voff_x, voff_c, voff_z = 0, n, 2 * n
# Offset of constraints.
numcon = 2 * n + 1
coff_bud, coff_abs1, coff_abs2 = 0, 1, 1 + n
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Variables (vector of x, c, z)
task.appendvars(numvar)
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (j + 1))
task.putvarname(voff_c + j, "c[%d]" % (j + 1))
task.putvarname(voff_z + j, "z[%d]" % (j + 1))
# Apply variable bounds (x >= 0, c and z free)
task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
task.putvarboundsliceconst(voff_c, voff_c + n, mosek.boundkey.fr, -inf, inf)
task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)
# Linear constraints
# - Budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
task.putaijlist([coff_bud] * n, range(voff_c, voff_c + n), m) # m^T c
rtemp = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp) # equals w + sum(x0)
# - Absolute value
task.appendcons(2 * n)
for i in range(0, n):
task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), n * [-1.0])
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), n * [1.0])
task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), n * [1.0])
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), n * [1.0])
task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)
# ACCs
aoff_q = 0
aoff_pow = k + 1
# - (gamma, GTx) in Q(k+1)
# The part of F and g for variable x:
# [0, 0, 0] [gamma]
# F = [GT, 0, 0], g = [0 ]
task.appendafes(k + 1)
task.putafeg(aoff_q, gamma)
for i in range(0, k):
task.putafefrow(aoff_q + 1 + i, range(voff_x, voff_x + n), GT[i])
qdom = task.appendquadraticconedomain(k + 1)
task.appendaccseq(qdom, aoff_q, None)
task.putaccname(aoff_q, "risk")
# - (c_j, 1, z_j) in P3(2/3, 1/3)
# The part of F and g for variables [c, z]:
# [0, I, 0] [0]
# F = [0, 0, I], g = [0]
# [0, 0, 0] [1]
task.appendafes(2 * n + 1)
task.putafefentrylist(range(aoff_pow, aoff_pow + n), range(voff_c, voff_c + n), [1.0] * n)
task.putafefentrylist(range(aoff_pow + n, aoff_pow + 2 * n), range(voff_z, voff_z + n), [1.0] * n)
task.putafeg(aoff_pow + 2 * n, 1.0)
# We use one row from F and g for both c_j and z_j, and the last row of F and g for the constant 1.
# NOTE: Here we reuse the last AFE and the power cone n times, but we store them only once.
powdom = task.appendprimalpowerconedomain(3, [2, 1])
afe_list = [(aoff_pow + i, aoff_pow + 2 * n, aoff_pow + n + i) for i in range(0, n)]
flat_afe_list = [idx for sublist in afe_list for idx in sublist]
task.appendaccs([powdom] * n, flat_afe_list, None)
for i in range(0, n):
task.putaccname(i + 1, "market_impact[%d]" % i)
# Objective
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Turn all log output off.
# task.putintparam(mosek.iparam.log,0)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
task.optimize()
# Display the solution summary for quick inspection of results.
task.solutionsummary(mosek.streamtype.msg)
# Check if the interior point solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itr)
if (solsta != mosek.solsta.optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
expret = 0.0
x = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
for j in range(0, n):
expret += mu[j] * x[j]
print("\nExpected return %e for gamma %e\n" % (expret, gamma))
portfolio_4_transcost.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_4_transcost.py
#
# Purpose : Implements a basic portfolio optimization model
# with fixed setup costs and transaction costs
# as a mixed-integer problem.
##
import mosek
import sys
if __name__ == '__main__':
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 ]
]
k = len(GT)
f = n*[0.01]
g = n*[0.001]
gamma = 0.36
# This value has no significance.
inf = 0.0
# Offset of variables.
numvar = 3 * n
voff_x, voff_z, voff_y = 0, n, 2 * n
# Offset of constraints.
numcon = 3 * n + 1
coff_bud, coff_abs1, coff_abs2, coff_swi = 0, 1, 1 + n, 1 + 2 * n
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Variables (vector of x, z, y)
task.appendvars(numvar)
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (j + 1))
task.putvarname(voff_z + j, "z[%d]" % (j + 1))
task.putvarname(voff_y + j, "y[%d]" % (j + 1))
# Apply variable bounds (x >= 0, z free, y binary)
task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)
task.putvarboundsliceconst(voff_y, voff_y + n, mosek.boundkey.ra, 0.0, 1.0)
task.putvartypelist(range(voff_y, voff_y + n), [mosek.variabletype.type_int] * n)
# Linear constraints
# - Budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([0] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
task.putaijlist([0] * n, range(voff_z, voff_z + n), g) # g^T z
task.putaijlist([0] * n, range(voff_y, voff_y + n), f) # f^T y
U = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, U, U) # equals w + sum(x0)
# - Absolute value
task.appendcons(2 * n)
for i in range(0, n):
task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), [-1.0] * n)
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), [1.0] * n)
task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), [1.0] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), [1.0] * n)
task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)
# - Switch
task.appendcons(n)
for i in range(0, n):
task.putconname(coff_swi + i, "switch[%d]" % (1 + i))
task.putaijlist(range(coff_swi, coff_swi + n), range(voff_z, voff_z + n), [1.0] * n)
task.putaijlist(range(coff_swi, coff_swi + n), range(voff_y, voff_y + n), [-U] * n)
task.putconboundlist(range(coff_swi, coff_swi + n), [mosek.boundkey.up] * n, [-inf] * n, [0.0] * n)
# ACCs
aoff_q = 0
# - (gamma, GTx) in Q(k+1)
# The part of F and g for variable x:
# [0, 0, 0] [gamma]
# F = [GT, 0, 0], g = [0 ]
task.appendafes(k + 1)
task.putafeg(aoff_q, gamma)
for i in range(0, k):
task.putafefrow(aoff_q + i + 1, range(voff_x, voff_x + n), GT[i])
qdom = task.appendquadraticconedomain(k + 1)
task.appendaccseq(qdom, aoff_q, None)
task.putaccname(0, "risk")
# Objective
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Turn all log output off.
# task.putintparam(mosek.iparam.log,0)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
task.optimize()
# Display the solution summary for quick inspection of results.
task.solutionsummary(mosek.streamtype.msg)
# Check if the integer solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itg)
if (solsta != mosek.solsta.integer_optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
expret = 0.0
x = task.getxxslice(mosek.soltype.itg, voff_x, voff_x + n)
for j in range(0, n):
expret += mu[j] * x[j]
tcost = 0.0
z = task.getxxslice(mosek.soltype.itg, voff_z, voff_z + n)
y = task.getxxslice(mosek.soltype.itg, voff_y, voff_y + n)
for j in range(0, n):
tcost += g[j] * z[j] + f[j] * y[j]
print("\nExpected return %e for gamma %e. Transaction cost: %e\n" % (expret, gamma, tcost))
portfolio_5_card.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_5_card.py
#
# Description : Implements a basic portfolio optimization model
# with cardinality constraints on number of assets traded.
##
import mosek
import sys
import numpy as np
# This value has no significance.
inf = 0.0
def markowitz_with_card(n, k, x0, w, gamma, mu, GT, K):
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Offset of variables.
numvar = 3 * n
voff_x, voff_z, voff_y = 0, n, 2 * n
# Offset of constraints.
numcon = 3 * n + 2
coff_bud, coff_abs1, coff_abs2, coff_swi, coff_card = 0, 1, 1 + n, 1 + 2 * n, 1 + 3 * n
# Variables (vector of x, z, y)
task.appendvars(numvar)
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (j + 1))
task.putvarname(voff_z + j, "z[%d]" % (j + 1))
task.putvarname(voff_y + j, "y[%d]" % (j + 1))
# Apply variable bounds (x >= 0, z free, y binary)
task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)
task.putvarboundsliceconst(voff_y, voff_y + n, mosek.boundkey.ra, 0.0, 1.0)
task.putvartypelist(range(voff_y, voff_y + n), [mosek.variabletype.type_int] * n)
# Linear constraints
# - Budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
U = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, U, U) # = w + sum(x0)
# - Absolute value
task.appendcons(2 * n)
for i in range(0, n):
task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), [-1.0] * n)
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), [1.0] * n)
task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), [1.0] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), [1.0] * n)
task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)
# - Switch
task.appendcons(n)
for i in range(0, n):
task.putconname(coff_swi + i, "switch[%d]" % (1 + i))
task.putaijlist(range(coff_swi, coff_swi + n), range(voff_z, voff_z + n), [1.0] * n)
task.putaijlist(range(coff_swi, coff_swi + n), range(voff_y, voff_y + n), [-U] * n)
task.putconboundlist(range(coff_swi, coff_swi + n), [mosek.boundkey.up] * n, [-inf] * n, [0.0] * n)
# - Cardinality
task.appendcons(1)
task.putconname(coff_card, "cardinality")
task.putaijlist([coff_card] * n, range(voff_y, voff_y + n), [1.0] * n) # e^T y
task.putconbound(coff_card, mosek.boundkey.up, -inf, K) # <= K
# ACCs
aoff_q = 0
# - (gamma, GTx) in Q(k+1)
# The part of F and g for variable x:
# [0, 0, 0] [gamma]
# F = [GT, 0, 0], g = [0 ]
task.appendafes(k + 1)
task.putafeg(aoff_q, gamma)
for i in range(0, k):
task.putafefrow(aoff_q + i + 1, range(voff_x, voff_x + n), GT[i])
qdom = task.appendquadraticconedomain(k + 1)
task.appendaccseq(qdom, aoff_q, None)
task.putaccname(0, "risk")
# Objective
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Turn all log output off.
task.putintparam(mosek.iparam.log,0)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
task.optimize()
# Check if the integer solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itg)
if (solsta != mosek.solsta.integer_optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
# Display the solution summary for quick inspection of results.
#task.solutionsummary(mosek.streamtype.msg)
return task.getxxslice(mosek.soltype.itg, voff_x + 0, voff_x + n)
if __name__ == '__main__':
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 ]
]
k = len(GT)
gamma = 0.25
for K in range(1, n+1):
xx = markowitz_with_card(n, k, x0, w, gamma, mu, GT, K)
expret = sum([xx[i]*mu[i] for i in range(n)])
print("Bound {0}: Return = {1:.4f} x = {2}".format(K, expret, xx))
portfolio_6_factor.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : portfolio_6_factor.py
#
# Purpose : Implements a basic portfolio optimization model
# of factor type.
##
import mosek
import sys
import numpy as np
if __name__ == '__main__':
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]
B = np.array([
[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 = np.array([
[0.0620, 0.0577],
[0.0577, 0.0908]
])
theta = np.array([0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459])
P = np.linalg.cholesky(S_F)
G_factor = B @ P
k = G_factor.shape[1]
gammas = [0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48]
inf = 0.0 # This value has no significance
# Variable offsets
numvar = n
voff_x = 0
# Constraints offsets
numcon = 1
coff_bud = 0
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Holding variable x of length n
# No other auxiliary variables are needed in this formulation
task.appendvars(numvar)
# Optionally we can give the variables names
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (1 + j))
# No short-selling in this model, all of x >= 0
task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)
# One linear constraint: total budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
rtemp = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp) # equals w + sum(x0)
# Input (gamma, G_factor_T x, diag(sqrt(theta))*x) in the AFE (affine expression) storage
# We need k+n+1 rows and we fill them in in three parts
task.appendafes(k + n + 1)
# 1. The first affine expression = gamma, will be specified later
# 2. The next k expressions comprise G_factor_T*x, we add them row by row
# transposing the matrix G_factor on the fly
for i in range(0, k):
task.putafefrow(i + 1, range(voff_x, voff_x + n), np.array(G_factor[:,i]))
# 3. The remaining n rows contain sqrt(theta) on the diagonal
task.putafefentrylist(range(k + 1, k + 1 + n), range(voff_x, voff_x + n), np.sqrt(theta))
# Input the affine conic constraint (gamma, G_factor_T x, diag(sqrt(theta))*x) \in QCone
# Add the quadratic domain of dimension k+n+1
qdom = task.appendquadraticconedomain(k + n + 1)
# Add the constraint
task.appendaccseq(qdom, 0, None)
task.putaccname(0, "risk")
# Objective: maximize expected return mu^T x
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
for gamma in gammas:
# Specify gamma in ACC
task.putafeg(0, gamma)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
# Solve the problem
task.optimize()
# Display solution summary for quick inspection of results.
# In this simplified example we skip checks for problem and solution status
task.solutionsummary(mosek.streamtype.msg)
# Check if the interior point solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itr)
if (solsta != mosek.solsta.optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
# Retrieve results
xx = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
expret = task.getprimalobj(mosek.soltype.itr)
print(f'Expected return: {expret:.10e} Std. deviation: {gamma:.4e}')
np.set_printoptions(precision=4)
print(f'Optimal portfolio: {np.array(xx)}')
pow1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : pow1.py
#
# 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
##
import sys
import mosek
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Only a symbolic constant
inf = 0.0
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
csub = [3, 4, 0]
cval = [1.0, 1.0, -1.0]
asub = [0, 1, 2]
aval = [1.0, 1.0, 0.5]
numvar, numcon = 5, 1 # x,y,z and 2 auxiliary variables for conic constraints
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
# Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
# Set up the linear part of the problem
task.putclist(csub, cval)
task.putarow(0, asub, aval)
task.putvarboundslice(0, numvar, [mosek.boundkey.fr] * numvar, [inf] * numvar, [inf] * numvar)
task.putconbound(0, mosek.boundkey.fx, 2.0, 2.0)
# Input the conic constraints
# Append two power cone domains
pc1 = task.appendprimalpowerconedomain(3, [0.2, 0.8])
pc2 = task.appendprimalpowerconedomain(3, [4.0, 6.0])
# Create data structures F,g so that
#
# F * x + g = (x(0), x(1), x(3), x(2), 1.0, x(4))
#
task.appendafes(6)
task.putafefentrylist([0, 1, 2, 3, 5], # Rows
[0, 1, 3, 2, 4], #Columns
[1.0, 1.0, 1.0, 1.0, 1.0])
task.putafeg(4, 1.0)
# Append the two conic constraints
task.appendacc(pc1, # Domain
[0, 1, 2], # Rows from F
None) # Unused
task.appendacc(pc2, # Domain
[3, 4, 5], # Rows from F
None) # Unused
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.maximize)
# Optimize the task
task.optimize()
# Print a summary containing information
# about the solution for debugging purposes
task.solutionsummary(mosek.streamtype.msg)
prosta = task.getprosta(mosek.soltype.itr)
solsta = task.getsolsta(mosek.soltype.itr)
# Output a solution
xx = task.getxx(mosek.soltype.itr)
if solsta == mosek.solsta.optimal:
print("Optimal solution: %s" % xx[0:3])
elif solsta == mosek.solsta.dual_infeas_cer:
print("Primal or dual infeasibility.\n")
elif solsta == mosek.solsta.prim_infeas_cer:
print("Primal or dual infeasibility.\n")
elif mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.code))
if msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
qcqo1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : qcqo1.py
#
# Purpose : Demonstrates how to solve small linear
# optimization problem using the MOSEK Python API.
##
import sys
import mosek
# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
# Set up and input bounds and linear coefficients
bkc = [mosek.boundkey.lo]
blc = [1.0]
buc = [inf]
bkx = [mosek.boundkey.lo,
mosek.boundkey.lo,
mosek.boundkey.lo]
blx = [0.0, 0.0, 0.0]
bux = [inf, inf, inf]
c = [0.0, -1.0, 0.0]
asub = [[0], [0], [0]]
aval = [[1.0], [1.0], [1.0]]
numvar = len(bkx)
numcon = len(bkc)
NUMANZ = 3
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
#Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
#Optionally add a constant term to the objective.
task.putcfix(0.0)
for j in range(numvar):
# Set the linear term c_j in the objective.
task.putcj(j, c[j])
# Set the bounds on variable j
# blx[j] <= x_j <= bux[j]
task.putvarbound(j, bkx[j], blx[j], bux[j])
# Input column j of A
task.putacol(j, # Variable (column) index.
# Row index of non-zeros in column j.
asub[j],
aval[j]) # Non-zero Values of column j.
for i in range(numcon):
task.putconbound(i, bkc[i], blc[i], buc[i])
# Set up and input quadratic objective
qsubi = [0, 1, 2, 2]
qsubj = [0, 1, 0, 2]
qval = [2.0, 0.2, -1.0, 2.0]
task.putqobj(qsubi, qsubj, qval)
# The lower triangular part of the Q^0
# matrix in the first constraint is specified.
# This corresponds to adding the term
# - x0^2 - x1^2 - 0.1 x2^2 + 0.2 x0 x2
qsubi = [0, 1, 2, 2]
qsubj = [0, 1, 2, 0]
qval = [-2.0, -2.0, -0.2, 0.2]
# put Q^0 in constraint with index 0.
task.putqconk(0, qsubi, qsubj, qval)
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.minimize)
# Optimize the task
task.optimize()
# Print a summary containing information
# about the solution for debugging purposes
task.solutionsummary(mosek.streamtype.msg)
prosta = task.getprosta(mosek.soltype.itr)
solsta = task.getsolsta(mosek.soltype.itr)
# Output a solution
xx = task.getxx(mosek.soltype.itr)
if solsta == mosek.solsta.optimal:
print("Optimal solution: %s" % xx)
elif solsta == mosek.solsta.dual_infeas_cer:
print("Primal or dual infeasibility.\n")
elif solsta == mosek.solsta.prim_infeas_cer:
print("Primal or dual infeasibility.\n")
elif mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
qo1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : qo1.py
#
# Purpose : Demonstrate how to solve a quadratic
# optimization problem using the MOSEK Python API.
##
import sys, os, mosek
# Since the actual value of Infinity is ignored, we define it solely
# for symbolic purposes:
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, streamprinter)
# Set up and input bounds and linear coefficients
bkc = [mosek.boundkey.lo]
blc = [1.0]
buc = [inf]
numvar = 3
bkx = [mosek.boundkey.lo] * numvar
blx = [0.0] * numvar
bux = [inf] * numvar
c = [0.0, -1.0, 0.0]
asub = [[0], [0], [0]]
aval = [[1.0], [1.0], [1.0]]
numvar = len(bkx)
numcon = len(bkc)
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
# Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
for j in range(numvar):
# Set the linear term c_j in the objective.
task.putcj(j, c[j])
# Set the bounds on variable j
# blx[j] <= x_j <= bux[j]
task.putvarbound(j, bkx[j], blx[j], bux[j])
# Input column j of A
task.putacol(j, # Variable (column) index.
# Row index of non-zeros in column j.
asub[j],
aval[j]) # Non-zero Values of column j.
for i in range(numcon):
task.putconbound(i, bkc[i], blc[i], buc[i])
# Set up and input quadratic objective
qsubi = [0, 1, 2, 2]
qsubj = [0, 1, 0, 2]
qval = [2.0, 0.2, -1.0, 2.0]
task.putqobj(qsubi, qsubj, qval)
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.minimize)
# Optimize
task.optimize()
# Print a summary containing information
# about the solution for debugging purposes
task.solutionsummary(mosek.streamtype.msg)
prosta = task.getprosta(mosek.soltype.itr)
solsta = task.getsolsta(mosek.soltype.itr)
# Output a solution
xx = task.getxx(mosek.soltype.itr)
if solsta == mosek.solsta.optimal:
print("Optimal solution: %s" % xx)
elif solsta == mosek.solsta.dual_infeas_cer:
print("Primal or dual infeasibility.\n")
elif solsta == mosek.solsta.prim_infeas_cer:
print("Primal or dual infeasibility.\n")
elif mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
if e.msg is not None:
import traceback
traceback.print_exc()
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
reoptimization.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : reoptimization.py
#
# Purpose: Demonstrates how to solve a linear
# optimization problem using the MOSEK API
# and modify and re-optimize the problem.
import sys
import mosek
# Since the actual value of Infinity is ignored, we define it solely
# for symbolic purposes:
inf = 0.0
def main():
# Create a task
with mosek.Task() as task:
# Bound keys for constraints
bkc = [mosek.boundkey.up,
mosek.boundkey.up,
mosek.boundkey.up]
# Bound values for constraints
blc = [-inf, -inf, -inf]
buc = [100000.0, 50000.0, 60000.0]
# Bound keys for variables
bkx = [mosek.boundkey.lo,
mosek.boundkey.lo,
mosek.boundkey.lo]
# Bound values for variables
blx = [0.0, 0.0, 0.0]
bux = [+inf, +inf, +inf]
# Objective coefficients
csub = [0, 1, 2]
cval = [1.5, 2.5, 3.0]
# We input the A matrix column-wise
# asub contains row indexes
asub = [0, 1, 2,
0, 1, 2,
0, 1, 2]
# acof contains coefficients
acof = [2.0, 3.0, 2.0,
4.0, 2.0, 3.0,
3.0, 3.0, 2.0]
# aptrb and aptre contains the offsets into asub and acof where
# columns start and end respectively
aptrb = [0, 3, 6]
aptre = [3, 6, 9]
numvar = len(bkx)
numcon = len(bkc)
# Append the constraints
task.appendcons(numcon)
# Append the variables.
task.appendvars(numvar)
# Input objective
task.putcfix(0.0)
task.putclist(csub, cval)
# Put constraint bounds
task.putconboundslice(0, numcon, bkc, blc, buc)
# Put variable bounds
task.putvarboundslice(0, numvar, bkx, blx, bux)
# Input A non-zeros by columns
for j in range(numvar):
ptrb, ptre = aptrb[j], aptre[j]
task.putacol(j,
asub[ptrb:ptre],
acof[ptrb:ptre])
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.maximize)
# Optimize the task
task.optimize()
# Output a solution
xx = task.getsolutionslice(mosek.soltype.bas,
mosek.solitem.xx,
0, numvar)
print("xx = {}".format(xx))
################# Make a change to the A matrix #############
task.putaij(0, 0, 3.0)
task.optimize()
# Output a solution
xx = task.getsolutionslice(mosek.soltype.bas,
mosek.solitem.xx,
0, numvar)
print("xx = {}".format(xx))
################### Add a new variable ######################
task.appendvars(1)
numvar+=1
# Set bounds on new varaible
task.putvarbound(task.getnumvar() - 1,
mosek.boundkey.lo,
0,
+inf)
# Change objective
task.putcj(task.getnumvar() - 1, 1.0)
# Put new values in the A matrix
acolsub = [0, 2]
acolval = [4.0, 1.0]
task.putacol(task.getnumvar() - 1, # column index
acolsub,
acolval)
# Change optimizer to simplex free and reoptimize
task.putintparam(mosek.iparam.optimizer,
mosek.optimizertype.free_simplex)
task.optimize()
# Output a solution
xx = task.getsolutionslice(mosek.soltype.bas,
mosek.solitem.xx,
0, numvar)
print("xx = {}".format(xx))
############# Add a new constraint #######################
task.appendcons(1)
numcon+=1
# Set bounds on new constraint
task.putconbound(task.getnumcon() - 1,
mosek.boundkey.up, -inf, 30000)
# Put new values in the A matrix
arowsub = [0, 1, 2, 3]
arowval = [1.0, 2.0, 1.0, 1.0]
task.putarow(task.getnumcon() - 1, # row index
arowsub,
arowval)
task.optimize()
# Output a solution
xx = task.getsolutionslice(mosek.soltype.bas,
mosek.solitem.xx,
0, numvar)
print("xx = {}".format(xx))
############# Change constraint bounds #######################
newbkc = [mosek.boundkey.up] * numcon
newblc = [-inf] * numcon
newbuc = [ 80000, 40000, 50000, 22000 ]
task.putconboundslice(0, numcon, newbkc, newblc, newbuc)
task.optimize()
# Output a solution
xx = task.getsolutionslice(mosek.soltype.bas,
mosek.solitem.xx,
0, numvar)
print("xx = {}".format(xx))
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
response.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : response.py
#
# Purpose : This example demonstrates proper response handling
# for problems solved with the interior-point optimizers.
#
##
import mosek
import sys
# A log message
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main(args):
filename = args[0] if len(args) >= 1 else "../data/cqo1.mps"
try:
with mosek.Task() as task:
# (Optional) set a log stream
# task.set_Stream(mosek.streamtype.log, streamprinter)
# (Optional) uncomment to see what happens when solution status is unknown
# task.putintparam(mosek.iparam.intpnt_max_iterations, 1)
# In this example we read data from a file
task.readdata(filename)
# Optimize
trmcode = task.optimize()
task.solutionsummary(mosek.streamtype.log)
# We expect solution status OPTIMAL
solsta = task.getsolsta(mosek.soltype.itr)
if solsta == mosek.solsta.optimal:
# Optimal solution. Fetch and print it.
print("An optimal interior-point solution is located.")
numvar = task.getnumvar()
xx = task.getxx(mosek.soltype.itr)
for i in range(numvar):
print("x[{0}] = {1}".format(i, xx[i]))
elif solsta == mosek.solsta.dual_infeas_cer:
print("Dual infeasibility certificate found.")
elif solsta == mosek.solsta.prim_infeas_cer:
print("Primal infeasibility certificate found.")
elif solsta == mosek.solsta.unknown:
# The solutions status is unknown. The termination code
# indicates why the optimizer terminated prematurely.
print("The solution status is unknown.")
symname, desc = mosek.Env.getcodedesc(trmcode)
print(" Termination code: {0} {1}".format(symname, desc))
else:
print("An unexpected solution status {0} is obtained.".format(str(solsta)))
except mosek.Error as e:
print("Unexpected error ({0}) {1}".format(e.errno, e.msg))
if __name__ == '__main__':
main(sys.argv[1:])
sdo1.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : sdo1.py
#
# Purpose : Demonstrates how to solve a small mixed semidefinite and conic quadratic
# optimization problem using the MOSEK Python API.
##
import sys
import mosek
# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task object and attach log stream printer
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, streamprinter)
# Bound keys for constraints
bkc = [mosek.boundkey.fx,
mosek.boundkey.fx]
# Bound values for constraints
blc = [1.0, 0.5]
buc = [1.0, 0.5]
# Below is the sparse representation of the A
# matrix stored by row.
asub = [[0],
[1, 2]]
aval = [[1.0],
[1.0, 1.0]]
barci = [0, 1, 1, 2, 2]
barcj = [0, 0, 1, 1, 2]
barcval = [2.0, 1.0, 2.0, 1.0, 2.0]
barai = [[0, 1, 2],
[0, 1, 2, 1, 2, 2]]
baraj = [[0, 1, 2],
[0, 0, 0, 1, 1, 2]]
baraval = [[1.0, 1.0, 1.0],
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]]
numvar = 3
numcon = len(bkc)
BARVARDIM = [3]
# Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
# Append 'numcon' empty constraints.
# The constraints will initially have no bounds.
task.appendcons(numcon)
# Append matrix variables of sizes in 'BARVARDIM'.
# The variables will initially be fixed at zero.
task.appendbarvars(BARVARDIM)
# Set the linear term c_0 in the objective.
task.putcj(0, 1.0)
for j in range(numvar):
# Set the bounds on variable j
# blx[j] <= x_j <= bux[j]
task.putvarbound(j, mosek.boundkey.fr, -inf, +inf)
for i in range(numcon):
# Set the bounds on constraints.
# blc[i] <= constraint_i <= buc[i]
task.putconbound(i, bkc[i], blc[i], buc[i])
# Input row i of A
task.putarow(i, # Constraint (row) index.
asub[i], # Column index of non-zeros in constraint i.
aval[i]) # Non-zero values of row i.
# Add the quadratic cone constraint
task.appendafes(3)
# Diagonal F matrix
task.putafefentrylist(range(3), range(3), [1.0]*3)
task.appendaccseq(task.appendquadraticconedomain(3), 0, None)
symc = task.appendsparsesymmat(BARVARDIM[0],
barci,
barcj,
barcval)
syma0 = task.appendsparsesymmat(BARVARDIM[0],
barai[0],
baraj[0],
baraval[0])
syma1 = task.appendsparsesymmat(BARVARDIM[0],
barai[1],
baraj[1],
baraval[1])
task.putbarcj(0, [symc], [1.0])
task.putbaraij(0, 0, [syma0], [1.0])
task.putbaraij(1, 0, [syma1], [1.0])
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.minimize)
# Solve the problem and print summary
task.optimize()
task.solutionsummary(mosek.streamtype.msg)
# Get status information about the solution
prosta = task.getprosta(mosek.soltype.itr)
solsta = task.getsolsta(mosek.soltype.itr)
if (solsta == mosek.solsta.optimal):
xx = task.getxx(mosek.soltype.itr)
barx = task.getbarxj(mosek.soltype.itr, 0)
print("Optimal solution:\nx=%s\nbarx=%s" % (xx, barx))
elif (solsta == mosek.solsta.dual_infeas_cer or
solsta == mosek.solsta.prim_infeas_cer):
print("Primal or dual infeasibility certificate found.\n")
elif solsta == mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
if e.msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
sdo2.py
#
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : sdo2.py
#
# 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.
from mosek import *
import sys
# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
# Sample data in sparse lower-triangular triplet form
C1_k = [0, 2]
C1_l = [0, 2]
C1_v = [1, 6]
A1_k = [0, 2, 2]
A1_l = [0, 0, 2]
A1_v = [1, 1, 2]
C2_k = [0, 1, 1, 2]
C2_l = [0, 0, 1, 2]
C2_v = [1, -3, 2, 1]
A2_k = [1, 1, 3]
A2_l = [0, 1, 3]
A2_v = [1, -1, -3]
b = 23
k = -3
# Create a task object and attach log stream printer
with Task() as task:
# Set log handler for debugging ootput
task.set_Stream(streamtype.log, streamprinter)
# Append two symmetric variables of dimension 3, 4
barvardims = [3, 4]
task.appendbarvars(barvardims)
# Semidefinite part of objective function
task.putbarcblocktriplet(
[0]*len(C1_v) + [1]*len(C2_v), # Which SDP variable (j)
C1_k + C2_k, # Entries: (k,l)->v
C1_l + C2_l,
C1_v + C2_v,
)
# Append two constraints
task.appendcons(2)
# First constraint (equality)
task.putbarablocktriplet(
[0]*(len(A1_v)+len(A2_v)), # Which constraint (i = 0)
[0]*len(A1_v) + [1]*len(A2_v), # Which SDP variable (j)
A1_k + A2_k, # Entries: (k,l)->v
A1_l + A2_l,
A1_v + A2_v,
)
# Second constraint (X2)_{1,2} <= k
task.putbarablocktriplet(
[1], # Which constraint (i = 1)
[1], # Which SDP variable (j = 1)
[1], [0], [0.5] # Entries: (k,l)->v
)
# Set bounds for constraints
task.putconboundlist([0,1], [boundkey.fx, boundkey.up],
[b, -inf],
[b, k])
# Write the problem for human inspection
task.writedata("test.ptf")
# Optimize
task.optimize()
task.solutionsummary(streamtype.msg)
# Get status information about the solution
solsta = task.getsolsta(soltype.itr)
if solsta == solsta.optimal:
# Assuming the optimization succeeded read solution
print("Solution (lower-triangular part vectorized): ")
for i in range(2):
X = task.getbarxj(soltype.itr, i)
print("X{i} = {X}".format(i=i, X=X))
elif (solsta == solsta.dual_infeas_cer or
solsta == solsta.prim_infeas_cer):
print("Primal or dual infeasibility certificate found.\n")
elif solsta == solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
sdo_lmi.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : sdo_lmi.py
#
# 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
import sys
from numpy import sqrt
import mosek
# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task object and attach log stream printer
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, streamprinter)
# Below is the sparse triplet representation of the F matrix.
afeidx = [0, 0, 1, 2, 2, 3]
varidx = [0, 1, 1, 0, 1, 0]
f_val = [-1, -1, 3, sqrt(2), sqrt(2), 3]
g = [0, -1, 0, -1]
barcj = [0, 0]
barck = [0, 1]
barcl = [0, 1]
barcv = [1, 1]
barfi = [0,0]
barfj = [0,0]
barfk = [0,1]
barfl = [0,0]
barfv = [0,1]
numvar = 2
numafe = 4
BARVARDIM = [2]
# Append 'numvar' variables.
# The variables will initially be fixed at zero (x=0).
task.appendvars(numvar)
# Append 'numafe' empty affine expressions.
task.appendafes(numafe)
# Append matrix variables of sizes in 'BARVARDIM'.
# The variables will initially be fixed at zero.
task.appendbarvars(BARVARDIM)
# Set the linear terms in the objective.
task.putcj(0, 1.0)
task.putcj(1, 1.0)
task.putcfix(1.0)
task.putbarcblocktriplet(barcj, barck, barcl, barcv)
for j in range(numvar):
# Set the bounds on variable j
# blx[j] <= x_j <= bux[j]
task.putvarbound(j, mosek.boundkey.fr, -inf, +inf)
# Set up the F matrix of the problem
task.putafefentrylist(afeidx, varidx, f_val)
# Set up the g vector of the problem
task.putafegslice(0, numafe, g)
task.putafebarfblocktriplet(barfi, barfj, barfk, barfl, barfv)
# Append R+ domain and the corresponding ACC
task.appendacc(task.appendrplusdomain(1), [0], None)
# Append SVEC_PSD domain and the corresponding ACC
task.appendacc(task.appendsvecpsdconedomain(3), [1,2,3], None)
# Input the objective sense (minimize/maximize)
task.putobjsense(mosek.objsense.minimize)
# Solve the problem and print summary
task.optimize()
task.solutionsummary(mosek.streamtype.msg)
# Get status information about the solution
prosta = task.getprosta(mosek.soltype.itr)
solsta = task.getsolsta(mosek.soltype.itr)
if (solsta == mosek.solsta.optimal):
xx = task.getxx(mosek.soltype.itr)
barx = task.getbarxj(mosek.soltype.itr, 0)
print("Optimal solution:\nx=%s\nbarx=%s" % (xx, barx))
elif (solsta == mosek.solsta.dual_infeas_cer or
solsta == mosek.solsta.prim_infeas_cer):
print("Primal or dual infeasibility certificate found.\n")
elif solsta == mosek.solsta.unknown:
print("Unknown solution status")
else:
print("Other solution status")
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
if e.msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
sensitivity.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : sensitivity.py
#
# Purpose : To demonstrate how to perform sensitivity
# analysis from the API on a small problem:
#
# minimize
#
# obj: +1 x11 + 2 x12 + 5 x23 + 2 x24 + 1 x31 + 2 x33 + 1 x34
# st
# c1: + x11 + x12 <= 400
# c2: + x23 + x24 <= 1200
# c3: + x31 + x33 + x34 <= 1000
# c4: + x11 + x31 = 800
# c5: + x12 = 100
# c6: + x23 + x33 = 500
# c7: + x24 + x34 = 500
#
# The example uses basis type sensitivity analysis.
##
import sys
import mosek
# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0
# Define a stream printer to grab output from MOSEK
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
# Create a task
with mosek.Task() as task:
# Attach a printer to the task
task.set_Stream(mosek.streamtype.log, streamprinter)
# Set up data
bkc = [mosek.boundkey.up, mosek.boundkey.up,
mosek.boundkey.up, mosek.boundkey.fx,
mosek.boundkey.fx, mosek.boundkey.fx,
mosek.boundkey.fx]
blc = [-inf, -inf, -inf, 800., 100., 500., 500.]
buc = [400., 1200., 1000., 800., 100., 500., 500.]
bkx = [mosek.boundkey.lo, mosek.boundkey.lo,
mosek.boundkey.lo, mosek.boundkey.lo,
mosek.boundkey.lo, mosek.boundkey.lo,
mosek.boundkey.lo]
c = [1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 1.0]
blx = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
bux = [inf, inf, inf, inf, inf, inf, inf]
ptrb = [0, 2, 4, 6, 8, 10, 12]
ptre = [2, 4, 6, 8, 10, 12, 14]
sub = [0, 3, 0, 4, 1, 5, 1, 6, 2, 3, 2, 5, 2, 6]
val = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
numcon = len(bkc)
numvar = len(bkx)
numanz = len(val)
# Input linear data
task.inputdata(numcon, numvar,
c, 0.0,
ptrb, ptre, sub, val,
bkc, blc, buc,
bkx, blx, bux)
# Set objective sense
task.putobjsense(mosek.objsense.minimize)
# Optimize
task.optimize()
# Analyze upper bound on c1 and the equality constraint on c4
subi = [0, 3]
marki = [mosek.mark.up, mosek.mark.up]
# Analyze lower bound on the variables x12 and x31
subj = [1, 4]
markj = [mosek.mark.lo, mosek.mark.lo]
(leftpricei, rightpricei, leftrangei, rightrangei,
leftpricej, rightpricej, leftrangej, rightrangej) = task.primalsensitivity(subi,
marki,
subj,
markj)
print('Results from sensitivity analysis on bounds:')
print('\tleftprice | rightprice | leftrange | rightrange ')
print('For constraints:')
for i in range(2):
print('\t%10f %10f %10f %10f' % (leftpricei[i],
rightpricei[i],
leftrangei[i],
rightrangei[i]))
print('For variables:')
for i in range(2):
print('\t%10f %10f %10f %10f' % (leftpricej[i],
rightpricej[i],
leftrangej[i],
rightrangej[i]))
subc = [2, 5]
leftprice, rightprice, leftrange, rightrange = task.dualsensitivity(subc)
print('Results from sensitivity analysis on objective coefficients:')
for i in range(2):
print('\t%10f %10f %10f %10f' % (leftprice[i],
rightprice[i],
leftrange[i],
rightrange[i]))
return None
# call the main function
try:
main()
except mosek.MosekException as e:
print("ERROR: %s" % str(e.errno))
if e.msg is not None:
print("\t%s" % e.msg)
sys.exit(1)
except:
import traceback
traceback.print_exc()
sys.exit(1)
simple.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : simple.py
#
# Purpose : Demonstrates a very simple example using MOSEK by
# reading a problem file, solving the problem and
# writing the solution to a file.
##
import mosek
import sys
def streamprinter(msg):
sys.stdout.write(msg)
sys.stdout.flush()
if len(sys.argv) <= 1:
print("Missing argument, syntax is:")
print(" simple inputfile [ solutionfile ]")
else:
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, streamprinter)
# We assume that a problem file was given as the first command
# line argument (received in `argv')
task.readdata(sys.argv[1])
# Solve the problem
task.optimize()
# Print a summary of the solution
task.solutionsummary(mosek.streamtype.log)
# If an output file was specified, save problem to a file
if len(sys.argv) >= 3:
# If using OPF format, these parameters will specify what to include in output
task.putintparam(mosek.iparam.opf_write_solutions, mosek.onoffkey.on)
task.putintparam(mosek.iparam.opf_write_problem, mosek.onoffkey.on)
task.putintparam(mosek.iparam.opf_write_hints, mosek.onoffkey.off)
task.putintparam(mosek.iparam.opf_write_parameters, mosek.onoffkey.off)
task.writedata(sys.argv[2])
solutionquality.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : solutionquality.py
#
# Purpose : To demonstrate how to examine the quality of a solution.
##
import sys
import mosek
def streamprinter(msg):
sys.stdout.write(msg)
sys.stdout.flush()
if len(sys.argv) <= 1:
print("Missing argument, syntax is:")
print(" solutionquality inputfile")
else:
try:
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, streamprinter)
# We assume that a problem file was given as the first command
# line argument (received in `argv')
task.readdata(sys.argv[1])
# Solve the problem
task.optimize()
# Print a summary of the solution
task.solutionsummary(mosek.streamtype.log)
whichsol = mosek.soltype.bas
solsta = task.getsolsta(whichsol)
pobj, pviolcon, pviolvar, pviolbarvar, pviolcones, pviolitg, \
dobj, dviolcon, dviolvar, dviolbarvar, dviolcones = \
task.getsolutioninfo(whichsol)
if solsta in [mosek.solsta.optimal]:
abs_obj_gap = abs(dobj - pobj)
rel_obj_gap = abs_obj_gap / \
(1.0 + min(abs(pobj), abs(dobj)))
max_primal_viol = max(pviolcon, pviolvar)
max_primal_viol = max(max_primal_viol, pviolbarvar)
max_primal_viol = max(max_primal_viol, pviolcones)
max_dual_viol = max(dviolcon, dviolvar)
max_dual_viol = max(max_dual_viol, dviolbarvar)
max_dual_viol = max(max_dual_viol, dviolcones)
# Assume the application needs the solution to be within
# 1e-6 ofoptimality in an absolute sense. Another approach
# would be looking at the relative objective gap
print("\n\n")
print("Customized solution information.\n")
print(" Absolute objective gap: %e\n" % abs_obj_gap)
print(" Relative objective gap: %e\n" % rel_obj_gap)
print(" Max primal violation : %e\n" % max_primal_viol)
print(" Max dual violation : %e\n" % max_dual_viol)
accepted = True
if rel_obj_gap > 1e-6:
print("Warning: The relative objective gap is LARGE.")
accepted = False
# We will accept a primal infeasibility of 1e-8 and
# dual infeasibility of 1e-6. These number should chosen problem
# dependent.
if max_primal_viol > 1e-8:
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:
numvar = task.getnumvar()
print("Optimal primal solution")
x = task.getxx(whichsol)
for j in range(numvar):
print("x[%d]: %e\n" % (j, x[j]))
else:
#Print detailed information about the solution
task.analyzesolution(mosek.streamtype.log, whichsol)
elif solsta in [mosek.solsta.dual_infeas_cer, mosek.solsta.prim_infeas_cer]:
print("Primal or dual infeasibility certificate found.")
elif solsta == mosek.solsta.unkwown:
print("The status of the solution is unknown.")
else:
print("Other solution status")
except mosek.Error as e:
print(e)
solvebasis.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : solvebasis.py
#
# Purpose : To demonstrate the usage of
# MSK_solvewithbasis on the problem:
#
# maximize x0 + x1
# st.
# x0 + 2.0 x1 <= 2
# x0 + x1 <= 6
# x0 >= 0, x1>= 0
#
# The problem has the slack variables
# xc0, xc1 on the constraints
# and the variabels x0 and x1.
#
# maximize x0 + x1
# st.
# x0 + 2.0 x1 -xc1 = 2
# x0 + x1 -xc2 = 6
# x0 >= 0, x1>= 0,
# xc1 <= 0 , xc2 <= 0
##
import mosek
import sys
def streamprinter(text):
sys.stdout.write(text)
sys.stdout.flush()
def main():
numcon = 2
numvar = 2
# Since the value infinity is never used, we define
# 'infinity' symbolic purposes only
infinity = 0
c = [1.0, 1.0]
ptrb = [0, 2]
ptre = [2, 3]
asub = [0, 1,
0, 1]
aval = [1.0, 1.0,
2.0, 1.0]
bkc = [mosek.boundkey.up,
mosek.boundkey.up]
blc = [-infinity,
-infinity]
buc = [2.0,
6.0]
bkx = [mosek.boundkey.lo,
mosek.boundkey.lo]
blx = [0.0,
0.0]
bux = [+infinity,
+infinity]
w1 = [2.0, 6.0]
w2 = [1.0, 0.0]
try:
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, streamprinter)
task.inputdata(numcon, numvar,
c,
0.0,
ptrb,
ptre,
asub,
aval,
bkc,
blc,
buc,
bkx,
blx,
bux)
task.putobjsense(mosek.objsense.maximize)
r = task.optimize()
if r != mosek.rescode.ok:
print("Mosek warning:", r)
basis = task.initbasissolve()
#List basis variables corresponding to columns of B
varsub = [0, 1]
for i in range(numcon):
if basis[varsub[i]] < numcon:
print("Basis variable no %d is xc%d" % (i, basis[i]))
else:
print("Basis variable no %d is x%d" %
(i, basis[i] - numcon))
# solve Bx = w1
# varsub contains index of non-zeros in b.
# On return b contains the solution x and
# varsub the index of the non-zeros in x.
nz = 2
nz = task.solvewithbasis(False, nz, varsub, w1)
print("nz = %s" % nz)
print("Solution to Bx = w1:")
for i in range(nz):
if basis[varsub[i]] < numcon:
print("xc %s = %s" % (basis[varsub[i]], w1[varsub[i]]))
else:
print("x%s = %s" %
(basis[varsub[i]] - numcon, w1[varsub[i]]))
# Solve B^Tx = w2
nz = 1
varsub[0] = 0
nz = task.solvewithbasis(True, nz, varsub, w2)
print("Solution to B^Tx = w2:")
for i in range(nz):
if basis[varsub[i]] < numcon:
print("xc %s = %s" % (basis[varsub[i]], w2[varsub[i]]))
else:
print("x %s = %s" %
(basis[varsub[i]] - numcon, w2[varsub[i]]))
except Exception as e:
print(e)
if __name__ == '__main__':
main()
solvelinear.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : solvelinear.py
#
# Purpose : To demonstrate the usage of MSK_solvewithbasis
# when solving the linear system:
#
# 1.0 x1 = b1
# -1.0 x0 + 1.0 x1 = b2
#
# with two different right hand sides
#
# b = (1.0, -2.0)
#
# and
#
# b = (7.0, 0.0)
##
import mosek
def setup(task,
aval,
asub,
ptrb,
ptre,
numvar,
basis):
# Since the value infinity is never used, we define
# 'infinity' symbolic purposes only
infinity = 0
skx = [mosek.stakey.bas] * numvar
skc = [mosek.stakey.fix] * numvar
task.appendvars(numvar)
task.appendcons(numvar)
for i in range(len(asub)):
task.putacol(i, asub[i], aval[i])
for i in range(numvar):
task.putconbound(i, mosek.boundkey.fx, 0.0, 0.0)
for i in range(numvar):
task.putvarbound(i,
mosek.boundkey.fr,
-infinity,
infinity)
# Define a basic solution by specifying
# status keys for variables & constraints.
task.deletesolution(mosek.soltype.bas);
task.putskcslice(mosek.soltype.bas, 0, numvar, skc);
task.putskxslice(mosek.soltype.bas, 0, numvar, skx);
task.initbasissolve(basis);
def main():
numcon = 2
numvar = 2
aval = [[-1.0],
[1.0, 1.0]]
asub = [[1],
[0, 1]]
ptrb = [0, 1]
ptre = [1, 3]
with mosek.Task() as task:
# Directs the log task stream to the user specified
# method task_msg_obj.streamCB
task.set_Stream(mosek.streamtype.log,
lambda msg: sys.stdout.write(msg))
# Put A matrix and factor A.
# Call this function only once for a given task.
basis = [0] * numvar
b = [0.0, -2.0]
bsub = [0, 1]
setup(task,
aval,
asub,
ptrb,
ptre,
numvar,
basis)
# now solve rhs
b = [1, -2]
bsub = [0, 1]
nz = task.solvewithbasis(False, 2, bsub, b)
print("\nSolution to Bx = b:\n")
# Print solution and show correspondents
# to original variables in the problem
for i in range(nz):
if basis[bsub[i]] < numcon:
print("This should never happen")
else:
print("x%d = %d" % (basis[bsub[i]] - numcon, b[bsub[i]]))
b[0] = 7
bsub[0] = 0
nz = task.solvewithbasis(False, 1, bsub, b)
print("\nSolution to Bx = b:\n")
# Print solution and show correspondents
# to original variables in the problem
for i in range(nz):
if basis[bsub[i]] < numcon:
print("This should never happen")
else:
print("x%d = %d" % (basis[bsub[i]] - numcon, b[bsub[i]]))
if __name__ == "__main__":
try:
main()
except:
import traceback
traceback.print_exc()
sparsecholesky.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: sparsecholesky.py
#
# Purpose: Demonstrate the sparse Cholesky factorization.
##
from __future__ import print_function
import mosek
import numpy
# Prints out a Cholesky factor presented in sparse form
def printsparse(n, perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc):
print("P = ", perm)
print("diag(D) = ", numpy.array(diag))
print("L=")
l = [[0.0 for j in range(n)] for i in range(n)]
for j in range(n):
for i in range(lptrc[j], lptrc[j] + lnzc[j]):
l[lsubc[i]][j] = lvalc[i]
print(numpy.array(l))
# Create the mosek environment.
with mosek.Env() as env:
# The matrix from the manual
# Observe that anzc, aptrc, asubc and avalc only specify the lower
# triangular part.
n = 4
anzc = [4, 1, 1, 1]
asubc = [0, 1, 2, 3, 1, 2, 3]
aptrc = [0, 4, 5, 6]
avalc = [4.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
b = [13.0, 3.0, 4.0, 5.0]
try:
perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc = env.computesparsecholesky(
0, # Mosek chooses number of threads
1, # Use reordering heuristic
1.0e-14,# Singularity tolerance
anzc, aptrc, asubc, avalc)
printsparse(n, perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc)
x = [b[p] for p in perm] # Permuted b is stored as x.
# Compute inv(L)*x.
env.sparsetriangularsolvedense(mosek.transpose.no,
lnzc, lptrc, lsubc, lvalc, x)
# Compute inv(L^T)*x.
env.sparsetriangularsolvedense(mosek.transpose.yes,
lnzc, lptrc, lsubc, lvalc, x)
print("\nSolution Ax=b: x = ", numpy.array(
[x[j] for i in range(n) for j in range(n) if perm[j] == i]), "\n")
except:
raise
#Example 2 - singular A
n = 3
anzc = [3, 2, 1]
asubc = [0, 1, 2, 1, 2, 2]
aptrc = [0, 3, 5]
avalc = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
try:
perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc = env.computesparsecholesky(
0, # Mosek chooses number of threads
1, # Use reordering heuristic
1.0e-14,# Singularity tolerance
anzc, aptrc, asubc, avalc)
printsparse(n, perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc)
except:
raise
