16 List of examples¶
List of examples shipped in the distribution of Fusion API for Python:
File |
Description |
|---|---|
Demonstrates a traffic network problem as a conic quadratic problem (CQO) |
|
A portfolio choice model |
|
A small bakery revenue maximization linear problem |
|
Shows how to break a long-running task |
|
An example of data/progress callback |
|
A simple conic exponential problem |
|
A simple conic quadratic problem |
|
Solving Stigler’s Nutrition model |
|
A simple problem with disjunctive constraints (DJC) |
|
Shows how to access the dual solution |
|
Linear regression with elastic net. Demonstrates model parametrization. |
|
Demonstrates a small one-facility location problem (CQO) |
|
A simple geometric program (GP) in conic form |
|
A simple linear problem |
|
Implements logistic regression and simple log-sum-exp (CEO) |
|
Computes the Lowner-John inner and outer ellipsoidal approximations of a polytope (SDO, Power Cone) |
|
Demonstrates how to solve the multi-processor scheduling problem and input an integer feasible point (MIP) |
|
A simple mixed-integer conic problem |
|
A simple mixed-integer linear problem |
|
A simple mixed-integer linear problem with an initial guess |
|
Solves the nearest correlation matrix problem (SDO, CQO) |
|
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 |
|
Implements a simple soft-margin Support Vector Machine (CQO) |
|
Demonstrate how to use SDP to solve convex relaxation of a mixed-integer QCQO 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 many matrix variables of the same dimension |
|
Models the cone of nonnegative polynomials and nonnegative trigonometric polynomials using Nesterov’s framework |
|
A SUDOKU solver (MIP) |
|
Demonstrates how to solve a total variation problem (CQO) |
|
Solves a simple Travelling Salesman Problem and shows how to add constraints to a model and re-optimize (MIP) |
Additional examples can be found on the MOSEK website and in other MOSEK publications.
TrafficNetworkModel.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: TrafficNetworkModel.py
#
# Purpose: Demonstrates a traffic network problem as a conic quadratic problem.
#
# Source: Robert Fourer, "Convexity Checking in Large-Scale Optimization",
# OR 53 --- Nottingham 6-8 September 2011.
#
# The problem:
# Given a directed graph representing a traffic network
# with one source and one sink, we have for each arc an
# associated capacity, base travel time and a
# sensitivity. Travel time along a specific arc increases
# as the flow approaches the capacity.
#
# Given a fixed inflow we now wish to find the
# configuration that minimizes the average travel time.
##
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
import sys
def main(args):
n = 4
arc_i = [0, 0, 2, 1, 2]
arc_j = [1, 2, 1, 3, 3]
arc_base = [4.0, 1.0, 2.0, 1.0, 6.0]
arc_cap = [10.0, 12.0, 20.0, 15.0, 10.0]
arc_sens = [0.1, 0.7, 0.9, 0.5, 0.1]
T = 20.0
source_idx = 0
sink_idx = 3
with Model() as M:
narcs = len(arc_i)
NxN = Set.make(n, n)
sens = Matrix.sparse(n, n, arc_i, arc_j, arc_sens)
cap = Matrix.sparse(n, n, arc_i, arc_j, arc_cap)
basetime = Matrix.sparse(n, n, arc_i, arc_j, arc_base)
cs_inv_matrix = \
Matrix.sparse(n, n, arc_i, arc_j,
[1.0 / (arc_sens[i] * arc_cap[i]) for i in range(narcs)])
s_inv_matrix = \
Matrix.sparse(n, n, arc_i, arc_j,
[1.0 / arc_sens[i] for i in range(narcs)])
sparsity = [ [ i, j] for i,j in zip(arc_i,arc_j) ]
x = M.variable("traffic_flow", [n,n], Domain.greaterThan(0.0).sparse(sparsity+[ [sink_idx, source_idx] ]))
t = M.variable("travel_time", [n,n], Domain.greaterThan(0.0).sparse(sparsity))
d = M.variable("d", [n,n], Domain.greaterThan(0.0).sparse(sparsity))
z = M.variable("z", [n,n], Domain.greaterThan(0.0).sparse(sparsity))
# Set the objective:
M.objective("Average travel time",
ObjectiveSense.Minimize,
(Expr.dot(basetime, x) + Expr.sum(d)) / T)
# Set up constraints
# Constraint (1a)
numnz = len(arc_sens)
v = Var.hstack([ d.pick(sparsity),
z.pick(sparsity),
x.pick(sparsity) ])
M.constraint("(1a)", v, Domain.inRotatedQCone(narcs,3))
# Constraint (1b)
M.constraint("(1b)", z + Expr.mulElm(x, cs_inv_matrix) == s_inv_matrix)
# Constraint (2)
M.constraint("(2)", Expr.sum(x,0) == Expr.sum(x,1))
# Constraint (3)
M.constraint("(3)", x[sink_idx, source_idx] == T)
M.writeTask("TrafficNetwork.ptf")
M.setLogHandler(sys.stdout);
M.solve()
M.writeTask("tn.ptf")
M.writeTask("tn.ptf")
flow = x.pick(sparsity).level()
print("Optimal flow:")
for i, j, v in zip(arc_i, arc_j, flow):
print("\tflow node%d->node%d = %f" % (i, j, v))
main(sys.argv[1:])
alan.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: alan.py
#
# Purpose: This file contains an implementation of the alan.gms (as
# found in the GAMS online model collection) using Fusion.
#
# The model is a simple portfolio choice model. The objective is to
# invest in a number of assets such that we minimize the risk, while
# requiring a certain expected return.
#
# We operate with 4 assets (hardware,software, show-biz and treasure
# bill). The risk is defined by the covariance matrix
# Q = [[ 4.0, 3.0, -1.0, 0.0 ],
# [ 3.0, 6.0, 1.0, 0.0 ],
# [ -1.0, 1.0, 10.0, 0.0 ],
# [ 0.0, 0.0, 0.0, 0.0 ]]
#
#
# We use the form Q = U^T * U, where U is a Cholesky factor of Q.
##
import sys
from mosek.fusion import *
import mosek.fusion.pythonic
#####################################################
### Problem data:
# Security names
securities = ["hardware", "software", "show-biz", "t-bills"]
# Two examples of mean returns on securities
mean1 = [9.0, 7.0, 11.0, 5.0]
mean2 = [8.0, 9.0, 12.0, 7.0]
# Target mean return
target = 10.0
# Factor of covariance matrix.
U = Matrix.dense([[2.0, 1.5, -0.5, 0.0],
[0.0, 1.93649167, 0.90369611, 0.0],
[0.0, 0.0, 2.98886824, 0.0],
[0.0, 0.0, 0.0, 0.0]])
numsec = len(securities)
###################################################
# Create a model with expected returns as parameter
with Model('alan') as M:
x = M.variable("x", numsec, Domain.greaterThan(0.0).withNamesOnAxis(securities,0))
t = M.variable("t", 1, Domain.greaterThan(0.0))
M.objective("minvar", ObjectiveSense.Minimize, t)
# sum securities to 1.0
M.constraint("wealth", Expr.sum(x) == 1.0)
# define target expected return
mean = M.parameter(numsec)
M.constraint("dmean", Expr.dot(mean, x) >= target)
# Bound on risk
M.constraint("t > ||Ux||^2", Expr.vstack(0.5, t, U@x) == Domain.inRotatedQCone())
M.writeTask('alan.ptf')
##################################################
# Define a method that solves a single instance
def solve(meanVal):
print("Solve using mean={0}".format(meanVal))
mean.setValue(meanVal)
M.solve()
solx = x.level()
print("Solution: {0}".format(solx))
return solx
#################################################
# Solve two instances with different means
sol1 = solve(mean1)
sol2 = solve(mean2)
baker.py
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: baker.py
#
# Purpose: Demonstrates a small linear problem.
#
# Source: "Linaer Algebra" by Knut Sydsaeter and Bernt Oeksendal.
#
# The problem: A baker has 150 kg flour, 22 kg sugar, 25 kg butter and two
# recipes:
# 1) Cakes, requiring 3.0 kg flour, 1.0 kg sugar and 1.2 kg butter per dozen.
# 2) Breads, requiring 5.0 kg flour, 0.5 kg sugar and 0.5 kg butter per dozen.
# Let the revenue per dozen cakes be $4 and the revenue per dozen breads be $6.
#
# We now wish to compute the combination of cakes and breads that will optimize
# the total revenue.
from mosek.fusion import *
import mosek.fusion.pythonic
ingredientnames = ["Flour", "Sugar", "Butter"]
stock = [150.0, 22.0, 25.0]
recipe = [[3.0, 5.0],
[1.0, 0.5],
[1.2, 0.5]]
productnames = ["Cakes", "Breads"]
revenue = [4.0, 6.0]
def main(args):
with Model("Recipe") as M:
# "production" defines the amount of each product to bake.
production = M.variable("production",
Set.make(productnames),
Domain.greaterThan(0.0))
# The objective is to maximize the total revenue.
M.objective("revenue",
ObjectiveSense.Maximize,
production.T @ revenue)
# The production is constrained by stock:
M.constraint(recipe @ production <= stock)
# We set stdout as the loghandler
M.setLogHandler(sys.stdout)
# We solve and fetch the solution:
M.solve()
print("Solution:")
res = production.level()
for name, val in zip(productnames, res):
print(" Number of %s : %.1f" % (name, val))
print(" Revenue : %.2f" % (res[0] * revenue[0] + res[1] * revenue[1]))
if __name__ == '__main__':
import sys
main(sys.argv[1:])
sys.exit(0)
breaksolver.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: breaksolver.py
#
# Purpose: Show how to break a long-running task.
##
import sys
from mosek.fusion import *
import random
import threading
import time
def main():
timeout = 5
n = 200 # number of binary variables
m = n // 5 # number of constraints
p = n // 5 # Each constraint picks p variables and requires that exactly half of them are 1
R = random.Random(1234)
print("Build problem...")
with Model('SolveBinary') as M:
M.setLogHandler(sys.stdout)
x = M.variable("x", n, Domain.binary())
M.objective(ObjectiveSense.Minimize, Expr.sum(x))
M.setLogHandler(sys.stdout)
L = list(range(n))
for i in range(m):
R.shuffle(L)
M.constraint(Expr.sum(x.pick(L[:p])), Domain.equalsTo(p // 2))
print("Start thread...")
T = threading.Thread(target=M.solve)
T0 = time.time()
try:
T.start() # optimization now running in background
# Loop until we get a solution or you run out of patience and press
# Ctrl-C
while True:
if not T.is_alive():
print("Solver terminated before anything happened!")
break
elif time.time() - T0 > timeout:
print("Solver terminated due to timeout!")
M.breakSolver()
break
except KeyboardInterrupt:
print("Signalling the solver that it can give up now!")
M.breakSolver()
finally:
try:
T.join() # wait for the solver to return
except:
pass
if __name__ == '__main__':
main()
callback.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
# File: callback.py
#
# Purpose: To demonstrate how to use the progress
# callback.
#
#
# callback psim
# callback dsim
# callback intpnt
# 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
from mosek.fusion import *
from mosek import callbackcode, iinfitem, dinfitem, liinfitem
def makeUserCallback(model, maxtime):
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.begin_bi:
print("Basis identification started.")
elif caller == callbackcode.end_bi:
print("Basis identification finished.")
else:
pass
if opttime >= maxtime:
# mosek is spending too much time. Terminate it.
print("Too much time, terminating.")
return 1
return 0
return userCallback
def userProgresCallback(caller):
# Handle the caller code here
pass
def main(slvr):
if slvr not in ['psim', 'dsim', 'intpnt']:
return
# Big random linear optimziation problem
n = 150
m = 700
A = numpy.random.uniform(0.0, 10.0, (m, n))
c = numpy.random.uniform(0.0, 10.0, n)
b = numpy.random.uniform(0.0, 10.0, m)
with Model() as M:
x = M.variable(n, Domain.unbounded())
M.constraint(Expr.mul(A, x), Domain.lessThan(b))
M.objective(ObjectiveSense.Maximize, Expr.dot(c, x))
if slvr == 'psim':
M.setSolverParam('optimizer', 'primalSimplex')
elif slvr == 'dsim':
M.setSolverParam('optimizer', 'dualSimplex')
else:
M.setSolverParam('optimizer', 'intpnt')
userCallback = makeUserCallback(model=M, maxtime=0.07)
M.setDataCallbackHandler(userCallback)
M.solve()
if __name__ == '__main__':
main(sys.argv[1] if len(sys.argv) > 1 else 'intpnt')
ceo1.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: ceo1.py
#
# Purpose: Demonstrates how to solve the problem
#
# minimize x1 + x2
# such that
# x1 + x2 + x3 = 1.0
# x1,x2 >= 0.0
# and x1 >= x2 * exp(x3/x2)
##
from mosek.fusion import *
import mosek.fusion.pythonic
with Model('ceo1') as M:
x = M.variable('x', 3, Domain.unbounded())
# Create the constraint
# x[0] + x[1] + x[2] = 1.0
M.constraint("lc", Expr.sum(x), Domain.equalsTo(1.0))
# Create the conic exponential constraint
expc = M.constraint("expc", x, Domain.inPExpCone())
# Set the objective function to (x[0] + x[1])
M.objective("obj", ObjectiveSense.Minimize, Expr.sum(x[0:2]))
# Solve the problem
M.solve()
M.writeTask("ceo1.ptf")
# Get the linear solution values
solx = x.level()
print('x1,x2,x3 = %s' % str(solx))
# Get conic solution of expc
expcval = expc.level()
expcdual = expc.dual()
print('expc levels = %s' % str(expcval))
print('expc dual conic var levels = %s' % str(expcdual))
cqo1.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: cqo1.py
#
# Purpose: Demonstrates how to solve the problem
#
# minimize y1 + y2 + y3
# such that
# x1 + x2 + 2.0 x3 = 1.0
# x1,x2,x3 >= 0.0
# and
# (y1,x1,x2) in C_3,
# (y2,y3,x3) in K_3
#
# where C_3 and K_3 are respectively the quadratic and
# rotated quadratic cone of size 3 defined as
# C_3 = { z1,z2,z3 : z1 >= sqrt(z2^2 + z3^2) }
# K_3 = { z1,z2,z3 : 2 z1 z2 >= z3^2 }
##
from mosek.fusion import *
import mosek.fusion.pythonic
with Model('cqo1') as M:
x = M.variable('x', 3, Domain.greaterThan(0.0))
y = M.variable('y', 3, Domain.unbounded())
# Create the aliases
# z1 = [ y[0],x[0],x[1] ]
# and z2 = [ y[1],y[2],x[2] ]
z1 = Var.vstack(y[0], x[0:2])
z2 = Var.vstack(y[1:3], x[2])
# Create the constraint
# x[0] + x[1] + 2.0 x[2] = 1.0
M.constraint("lc", Expr.dot([1.0, 1.0, 2.0], x) == 1.0)
# Create the constraints
# z1 belongs to C_3
# z2 belongs to K_3
# where C_3 and K_3 are respectively the quadratic and
# rotated quadratic cone of size 3, i.e.
# z1[0] >= sqrt(z1[1]^2 + z1[2]^2)
# and 2.0 z2[0] z2[1] >= z2[2]^2
qc1 = M.constraint("qc1", z1, Domain.inQCone())
qc2 = M.constraint("qc2", z2, Domain.inRotatedQCone())
# Set the objective function to (y[0] + y[1] + y[2])
M.objective("obj", ObjectiveSense.Minimize, Expr.sum(y))
# Solve the problem
M.solve()
M.writeTask('cqo1.ptf')
# Get the linear solution values
solx = x.level()
soly = y.level()
print('x1,x2,x3 = %s' % str(solx))
print('y1,y2,y3 = %s' % str(soly))
# Get conic solution of qc1
qc1lvl = qc1.level()
qc1sn = qc1.dual()
print('qc1 levels = %s' % str(qc1lvl))
print('qc1 dual conic var levels = %s' % str(qc1sn))
diet.py
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: diet.py
#
# Purpose: Solving Stigler's Nutrition model (DIET,SEQ=7)
#
# Source: GAMS Model library,
# Dantzig, G B, Chapter 27.1. In Linear Programming and Extensions.
# Princeton University Press, Princeton, New Jersey, 1963.
#
# Given a set of nutrients, foods with associated nutrient values, allowance of
# nutrients per day, the model will find the cheapest combination of foods
# which will provide the necessary nutrients.
from mosek.fusion import Model, Matrix, Domain, Expr, ObjectiveSense, Set
import numpy
##
# Arguments for construction:
# name - Model name.
# foods - List of M names of the foods.
# nutrients - List of N names of the nutrients.
# daily_allowance - List of N floats denoting the daily allowance of each
# nutrient.
# nutritive_value - Two-dimensional MxN array of floats where each row
# denotes the nutrient values for a single food per $ spent.
def DietModel(name,
foods,
nutrients,
daily_allowance,
nutritive_value):
with Model() as DM:
dailyAllowance = numpy.array(daily_allowance, float)
nutrientValue = Matrix.dense(nutritive_value).transpose()
M = len(foods)
N = len(nutrients)
if len(dailyAllowance) != N:
raise ValueError("Length of daily_allowance does not match "
"the number of nutrients")
if nutrientValue.numColumns() != M:
raise ValueError("Number of rows in nutrient_value does not "
"match the number of foods")
if nutrientValue.numRows() != N:
raise ValueError("Number of columns in nutrient_value does "
"not match the number of nutrients")
dailyPurchase = DM.variable(
'Daily Purchase',
[M],
Domain.greaterThan(0.0).withNamesOnAxis(foods,0))
dailyNutrients = DM.constraint(
'Nutrient Balance',
Expr.mul(nutrientValue,
dailyPurchase),
Domain.greaterThan(dailyAllowance).withNamesOnAxis(nutrients,0))
DM.objective(ObjectiveSense.Minimize,
Expr.sum(dailyPurchase))
DM.solve()
DM.writeTask("diet.ptf")
return dailyPurchase.level(),dailyNutrients.level()
# class DietModel(Model):
# def __init__(self, name,
# foods,
# nutrients,
# daily_allowance,
# nutritive_value):
# Model.__init__(self, name)
# finished = False
# try:
# self.foods = numpy.array([str(f) for f in foods],dtype=numpy.dtype(object))
# print(self.foods,len(numpy.array(self.foods,dtype=numpy.dtype(object))))
# self.nutrients = [str(n) for n in nutrients]
# self.dailyAllowance = numpy.array(daily_allowance, float)
# self.nutrientValue = Matrix.dense(nutritive_value).transpose()
# M = len(self.foods)
# N = len(self.nutrients)
# if len(self.dailyAllowance) != N:
# raise ValueError("Length of daily_allowance does not match "
# "the number of nutrients")
# if self.nutrientValue.numColumns() != M:
# raise ValueError("Number of rows in nutrient_value does not "
# "match the number of foods")
# if self.nutrientValue.numRows() != N:
# raise ValueError("Number of columns in nutrient_value does "
# "not match the number of nutrients")
# self.__dailyPurchase = self.variable('Daily Purchase',
# [M],
# Domain.greaterThan(0.0).withNamesOnAxis(self.foods,0))
# self.__dailyNutrients = \
# self.constraint('Nutrient Balance',
# Expr.mul(self.nutrientValue,
# self.__dailyPurchase),
# Domain.greaterThan(self.dailyAllowance).withNamesOnAxis(self.nutrients,0))
# self.objective(ObjectiveSense.Minimize,
# Expr.sum(self.__dailyPurchase))
# finished = True
# finally:
# if not finished:
# self.__del__()
# def dailyPurchase(self):
# return self.__dailyPurchase.level()
# def dailyNutrients(self):
# return self.__dailyNutrients.level()
if __name__ == '__main__':
nutrient_unit = [
"thousands", "grams", "grams",
"milligrams", "thousand ius", "milligrams",
"milligrams", "milligrams", "milligrams"]
nutrients = [
"calorie", "protein", "calcium",
"iron", "vitamin a", "vitamin b1",
"vitamin b2", "niacin", "vitamin c"]
foods = [
"wheat", "cornmeal", "cannedmilk", "margarine", "cheese",
"peanut butter", "lard", "liver", "porkroast", "salmon",
"greenbeans", "cabbage", "onions", "potatoes", "spinach",
"sweet potatos", "peaches", "prunes", "limabeans", "navybeans"]
nutritive_value = [
# calorie calcium vitamin a vitamin b2 vitamin c
# protein iron vitamin b1 niacin
[44.7, 1411, 2.0, 365, 0, 55.4, 33.3, 441, 0], # wheat
[36, 897, 1.7, 99, 30.9, 17.4, 7.9, 106, 0], # cornmeal
[8.4, 422, 15.1, 9, 26, 3, 23.5, 11, 60], # cannedmilk
[20.6, 17, .6, 6, 55.8, .2, 0, 0, 0], # margarine
[7.4, 448, 16.4, 19, 28.1, .8, 10.3, 4, 0], # cheese
[15.7, 661, 1, 48, 0, 9.6, 8.1, 471, 0], # peanut butter
[41.7, 0, 0, 0, .2, 0, .5, 5, 0], # lard
[2.2, 333, .2, 139, 169.2, 6.4, 50.8, 316, 525], # liver
[4.4, 249, .3, 37, 0, 18.2, 3.6, 79, 0], # porkroast
[5.8, 705, 6.8, 45, 3.5, 1, 4.9, 209, 0], # salmon
[2.4, 138, 3.7, 80, 69, 4.3, 5.8, 37, 862], # greenbeans
[2.6, 125, 4, 36, 7.2, 9, 4.5, 26, 5369], # cabbage
[5.8, 166, 3.8, 59, 16.6, 4.7, 5.9, 21, 1184], # onions
[14.3, 336, 1.8, 118, 6.7, 29.4, 7.1, 198, 2522], # potatoes
[1.1, 106, 0.0, 138, 918.4, 5.7, 13.8, 33, 2755], # spinach
[9.6, 138, 2.7, 54, 290.7, 8.4, 5.4, 83, 1912], # sweet potatos
[8.5, 87, 1.7, 173, 86.8, 1.2, 4.3, 55, 57], # peaches
[12.8, 99, 2.5, 154, 85.7, 3.9, 4.3, 65, 257], # prunes
[17.4, 1055, 3.7, 459, 5.1, 26.9, 38.2, 93, 0], # limabeans
[26.9, 1691, 11.4, 792, 0, 38.4, 24.6, 217, 0]] # navybeans
daily_allowance = [
3., 70., .8, 12., 5., 1.8, 2.7, 18., 75.]
x,y = DietModel('Stinglers Diet Problem',
foods,
nutrients,
daily_allowance,
nutritive_value)
# with DietModel('Stinglers Diet Problem',
# foods,
# nutrients,
# daily_allowance,
# nutritive_value) as M:
# M.solve()
if True:
print("Solution:")
#x = M.dailyPurchase()
for name, quantum in zip(foods, x):
print("%20s : $%.4f " % (name, quantum))
print("Nutrients:")
#y = M.dailyNutrients()
for name, quantum, unit, minval in zip(nutrients, y, nutrient_unit, daily_allowance):
print("%20s : %7.2f %s (%.2f)" % (name, quantum, unit, minval))
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
##
from mosek.fusion import *
import mosek.fusion.pythonic
import sys
with Model('djc1') as M:
# Create a variable of length 4
x = M.variable('x', 4)
# First disjunctive constraint
M.disjunction( (x[0]-2*x[1]<=-1) & (x[2:4]==0)
|
(x[2]-3*x[3]<=-2) & (x[0:2]==0))
# Second disjunctive constraint
# Array of terms reading x_i = 2.5 for i = 0,1,2,3
M.disjunction( [ x[i] == 2.5 for i in range(4) ])
# The linear constraint
M.constraint(Expr.sum(x) >= -10)
# Objective
M.objective(ObjectiveSense.Minimize, Expr.dot([2,1,3,1], x))
# Useful for debugging
M.writeTask("djc.ptf") # Save to a readable file
M.setLogHandler(sys.stdout) # Enable log output
# Solve
M.solve()
if M.getPrimalSolutionStatus() == SolutionStatus.Optimal:
print(f"Solution ={x.level()}")
else:
print("Another solution status")
duality.py
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: duality.py
#
# Purpose: Shows how to access the dual solution
import sys
import mosek
import mosek.fusion
from mosek.fusion import *
def main(args):
A = [[-0.5, 1.0]]
b = [1.0]
c = [1.0, 1.0]
with Model("duality1") as M:
x = M.variable("x", 2, Domain.greaterThan(0.0))
con = M.constraint(
Expr.sub(Expr.mul(Matrix.dense(A), x), b), Domain.equalsTo(0.0))
M.objective("obj", ObjectiveSense.Minimize, Expr.dot(c, x))
M.solve()
print("x =", x.level())
print("s =", x.dual())
print("y =", con.dual())
if __name__ == '__main__':
main(sys.argv[1:])
elastic.py
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: elastic.py
#
# Purpose: Demonstrates model parametrization on the example of an elastic net linear regression:
#
# min_x |Ax-b|_2 + lambda1*|x|_1 + lambda2*|x|_2
from mosek.fusion import *
import mosek.fusion.pythonic
# Construct the model with parameters b, lambda1, lambda2
# and with prescribed matrix A
def initializeModel(m, n, A):
M = Model()
x = M.variable("x", n)
# t >= |Ax-b|_2 where b is a parameter
b = M.parameter("b", m)
t = M.variable()
M.constraint(Expr.vstack(t, A@x-b) == Domain.inQCone())
# p_i >= |x_i|, i=1..n
p = M.variable(n)
M.constraint(Expr.hstack(p, x) == Domain.inQCone())
# q >= |x|_2
q = M.variable()
M.constraint(Expr.vstack(q, x) == Domain.inQCone())
# Objective, parametrized with lambda1, lambda2
# t + lambda1*sum(p) + lambda2*q
lambda1 = M.parameter("lambda1")
lambda2 = M.parameter("lambda2")
obj = t + lambda1 * Expr.sum(p) + lambda2 * q
M.objective(ObjectiveSense.Minimize, obj)
# Return the ready model
return M
def smallExample():
# Create a small example
m, n = 4, 2
A = [ [1.0, 2.0],
[3.0, 4.0],
[-2.0, -1.0],
[-4.0, -3.0] ]
M = initializeModel(m, n, A)
# For convenience retrieve some elements of the model
b = M.getParameter("b")
lambda1 = M.getParameter("lambda1")
lambda2 = M.getParameter("lambda2")
x = M.getVariable("x")
# First solve
b.setValue([0.1, 1.2, -1.1, 3.0])
lambda1.setValue(0.1)
lambda2.setValue(0.01)
M.solve()
print("Objective {0}, solution x = {1}".format(M.primalObjValue(), x.level()))
# Increase lambda1
lambda1.setValue(0.5)
M.solve()
print("Objective {0}, solution x = {1}".format(M.primalObjValue(), x.level()))
# Now change the data completely
b.setValue([1.0, 1.0, 1.0, 1.0])
lambda1.setValue(0.0)
lambda2.setValue(0.0)
M.solve()
print("Objective {0}, solution x = {1}".format(M.primalObjValue(), x.level()))
# And increase lamda2
lambda2.setValue(1.4145)
M.solve()
print("Objective {0}, solution x = {1}".format(M.primalObjValue(), x.level()))
M.dispose()
smallExample()
facility_location.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: facility_location.py
#
# Purpose: Demonstrates a small one-facility location problem.
#
# Given 10 customers placed in a grid we wish to place a facility
# somewhere so that the total sum of distances to customers is
# minimized.
#
# The problem is formulated as a conic optimization problem as follows.
# Let f=(fx,fy) be the (unknown) location of the facility, and let
# c_i=(cx_i,cy_i) be the (known) customer locations; then we wish to
# minimize
# sum_i || f - c_i ||
# where
# ||.||
# denotes the euclidian norm.
# This is formulated as
#
# minimize sum(d_i)
# such that d_i ^ 2 > tx_i ^ 2 + ty_i ^ 2, for all i (a.k.a. (d_i,tx_i,ty_i) in C^3_r)
# tx_i = cx_i - fx, for all i
# ty_i = cy_i - fy, for all i
# d_i > 0, for all i
##
import sys
from mosek.fusion import *
# Customer locations
customerloc = Matrix.dense([[12.0, 2.0],
[15.0, 13.0],
[10.0, 8.0],
[0.0, 10.0],
[6.0, 13.0],
[5.0, 8.0],
[10.0, 12.0],
[4.0, 6.0],
[5.0, 2.0],
[1.0, 10.0]])
N = customerloc.numRows()
with Model('FacilityLocation') as M:
# Variable holding the facility location
f = M.variable("facility", Set.make(1, 2), Domain.unbounded())
# Variable defining the euclidian distances to each customer
d = M.variable("dist", Set.make(N, 1), Domain.greaterThan(0.0))
# Variable defining the x and y differences to each customer
t = M.variable("t", Set.make(N, 2), Domain.unbounded())
M.constraint('dist measure', Var.hstack(d, t), Domain.inQCone())
fxy = Var.repeat(f, N)
M.constraint("xy diff", Expr.add(t, fxy), Domain.equalsTo(customerloc))
M.objective("total_dist", ObjectiveSense.Minimize, Expr.sum(d))
M.solve()
M.writeTask("facility_location.task")
print('Facility location =', f.level())
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, array
from mosek.fusion import *
import mosek.fusion.pythonic
import sys
# Models log(sum(exp(Ax+b))) <= 0.
# Each row of [A b] describes one of the exp-terms
def logsumexp(M, A, x, b):
k = int(A.shape[0])
u = M.variable(k)
M.constraint(Expr.sum(u) == 1.0)
M.constraint(Expr.hstack(u,
Expr.constTerm(k, 1.0),
A @ x + b), Domain.inPExpCone())
# 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):
with Model('max_vol_box') as M:
xyz = M.variable(3)
M.objective('Objective', ObjectiveSense.Maximize, Expr.sum(xyz))
logsumexp(M, array([[1,1,0],[1,0,1]]), xyz, array([log(2.0/Aw), log(2.0/Aw)]))
M.constraint(Expr.dot([0, 1, 1], xyz) <= log(Af))
M.constraint(Expr.dot([1,-1, 0], xyz) == Domain.inRange(log(alpha),log(beta)))
M.constraint(Expr.dot([0,-1, 1], xyz) == Domain.inRange(log(gamma),log(delta)))
M.setLogHandler(sys.stdout)
M.solve()
return exp(xyz.level())
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))
lo1.py
####
## Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
##
## File: lo1.py
##
## Purpose: Demonstrates how to solve the problem
##
## maximize 3*x0 + 1*x1 + 5*x2 + x3
## such that
## 3*x0 + 1*x1 + 2*x2 = 30,
## 2*x0 + 1*x1 + 3*x2 + 1*x3 > 15,
## 2*x1 + + 3*x3 < 25
## and
## x0,x1,x2,x3 > 0,
## 0 < x1 < 10
####
import sys
from mosek.fusion import *
import mosek.fusion.pythonic
def main(args):
A = [[3.0, 1.0, 2.0, 0.0],
[2.0, 1.0, 3.0, 1.0],
[0.0, 2.0, 0.0, 3.0]]
c = [3.0, 1.0, 5.0, 1.0]
# Create a model with the name 'lo1'
with Model("lo1") as M:
# Create variable 'x' of length 4
x = M.variable("x", 4, Domain.greaterThan(0.0))
# Create constraints
M.constraint(x[1] <= 10)
M.constraint("c1", x.T @ A[0] == 30.0)
M.constraint("c2", x.T @ A[1] >= 15.0)
M.constraint("c3", x.T @ A[2] <= 25.0)
# Set the objective function to (c^t * x)
M.objective("obj", ObjectiveSense.Maximize, x.T @ c)
# Solve the problem
M.solve()
# Get the solution values
sol = x.level()
print('\n'.join(["x[%d] = %f" % (i, sol[i]) for i in range(4)]))
if __name__ == '__main__':
main(sys.argv[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 Fusion.
#
# Plots an example for 2D datasets.
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
import sys, itertools
# t >= log( 1 + exp(u) ) coordinatewise
def softplus(M, t, u):
n = t.getShape()[0]
z1 = M.variable(n)
z2 = M.variable(n)
M.set(z1 + z2 == 1,
Expr.hstack(z1, Expr.constTerm(n, 1.0), u-t) == Domain.inPExpCone(),
Expr.hstack(z2, Expr.constTerm(n, 1.0), -t) == Domain.inPExpCone())
# 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
M = Model()
theta = M.variable(d)
t = M.variable(n)
reg = M.variable()
M.objective(ObjectiveSense.Minimize, Expr.sum(t) + lamb * reg)
M.constraint(Var.vstack(reg, theta), Domain.inQCone())
signs = list(map(lambda y: -1.0 if y==1 else 1.0, y))
softplus(M, t, Expr.mulElm(X @ theta, signs))
return M, theta
# 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)
M, theta = logisticRegression(X, y, lamb)
M.solve()
plot2d(x, y, 6, theta.level())
'''
# 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)
M, theta = logisticRegression(X, y, 0.1)
M.setLogHandler(sys.stdout)
M.solve()
try:
print(theta.level())
except:
print("Could not get solution")
sys.exit(1)
lownerjohn_ellipsoid.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: lownerjohn_ellipsoid.py
#
# Purpose:
# Computes the Lowner-John inner and outer ellipsoidal
# approximations of a polytope.
#
# Note:
# To plot the solution the Python package matplotlib is required.
#
# References:
# [1] "Lectures on Modern Optimization", Ben-Tal and Nemirovski, 2000.
# [2] "MOSEK modeling manual", 2018
##
import sys
from math import sqrt, ceil, log
from mosek.fusion import *
import mosek.fusion.pythonic
'''
Purpose: Models the hypograph of the n-th power of the
determinant of a positive definite matrix. See [1,2] for more details.
The convex set (a hypograph)
C = { (X, t) \\in S^n_+ x R | t <= det(X)^{1/n} },
can be modeled as the intersection of a semidefinite cone
[ X, Z; Z^T Diag(Z) ] >= 0
and a geometric mean bound
t <= (Z11*Z22*...*Znn)^{1/n}
'''
def det_rootn(M, t, n):
# Setup variables
Y = M.variable(Domain.inPSDCone(2 * n))
# Setup Y = [X, Z; Z^T , diag(Z)]
X = Y[0:n, 0:n]
Z = Y[0:n, n:2*n]
DZ = Y[n:2*n, n:2*n]
# Z is lower-triangular
M.constraint(Z.pick([[i,j] for i in range(n) for j in range(i+1,n)]), Domain.equalsTo(0.0))
# DZ = Diag(Z)
M.constraint(DZ == Expr.mulElm(Z, Matrix.eye(n)))
# (Z11*Z22*...*Znn) >= t^n
M.constraint(Expr.vstack(DZ.diag(), t), Domain.inPGeoMeanCone())
# Return an n x n PSD variable which satisfies t <= det(X)^(1/n)
return X
'''
The inner ellipsoidal approximation to a polytope
S = { x \\in R^n | Ax < b }.
maximizes the volume of the inscribed ellipsoid,
{ x | x = C*u + d, || u ||_2 <= 1 }.
The volume is proportional to det(C)^(1/n), so the
problem can be solved as
maximize t
subject to t <= det(C)^(1/n)
|| C*ai ||_2 <= bi - ai^T * d, i=1,...,m
C is PSD
which is equivalent to a mixed conic quadratic and semidefinite
programming problem.
'''
def lownerjohn_inner(A, b):
with Model("lownerjohn_inner") as M:
M.setLogHandler(sys.stdout)
m, n = len(A), len(A[0])
# Setup variables
t = M.variable("t", 1, Domain.greaterThan(0.0))
C = det_rootn(M, t, n)
d = M.variable("d", n, Domain.unbounded())
# (b-Ad, AC) generate cones
M.constraint("qc", Expr.hstack(b - A @ d, A @ C), Domain.inQCone())
# Objective: Maximize t
M.objective(ObjectiveSense.Maximize, t)
M.writeTask("lj-inner.ptf")
M.solve()
C, d = C.level(), d.level()
return ([C[i:i + n] for i in range(0, n * n, n)], d)
'''
The outer ellipsoidal approximation to a polytope given
as the convex hull of a set of points
S = conv{ x1, x2, ... , xm }
minimizes the volume of the enclosing ellipsoid,
{ x | || P*x-c ||_2 <= 1 }
The volume is proportional to det(P)^{-1/n}, so the problem can
be solved as
maximize t
subject to t <= det(P)^(1/n)
|| P*xi - c ||_2 <= 1, i=1,...,m
P is PSD.
'''
def lownerjohn_outer(x):
with Model("lownerjohn_outer") as M:
M.setLogHandler(sys.stdout)
m, n = len(x), len(x[0])
# Setup variables
t = M.variable("t", 1, Domain.greaterThan(0.0))
P = det_rootn(M, t, n)
c = M.variable("c", Domain.unbounded().withShape(1,n))
# (1, Px-c) in cone
M.constraint("qc", Expr.hstack(Expr.ones(m), x @ P - Expr.repeat(c,m,0)), Domain.inQCone())
# Objective: Maximize t
M.objective(ObjectiveSense.Maximize, t)
M.writeTask("lj-outer.ptf")
M.solve()
P, c = P.level(), c.level()
return ([P[i:i + n] for i in range(0, n * n, n)], c)
##########################################################################
if __name__ == '__main__':
#Vertices of a polygon in 2D
p = [[0., 0.], [1., 3.], [5.5, 4.5], [7., 4.], [7., 1.], [3., -2.]]
nVerts = len(p)
#The hyperplane representation of the same polytope
A = [[-p[i][1] + p[i - 1][1], p[i][0] - p[i - 1][0]]
for i in range(len(p))]
b = [A[i][0] * p[i][0] + A[i][1] * p[i][1] for i in range(len(p))]
Po, co = lownerjohn_outer(p)
Ci, di = lownerjohn_inner(A, b)
#Visualization
try:
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as patches
#Polygon
fig = plt.figure()
ax = fig.add_subplot(111)
ax.add_patch(patches.Polygon(p, fill=False, color="red"))
#The inner ellipse
theta = np.linspace(0, 2 * np.pi, 100)
x = Ci[0][0] * np.cos(theta) + Ci[0][1] * np.sin(theta) + di[0]
y = Ci[1][0] * np.cos(theta) + Ci[1][1] * np.sin(theta) + di[1]
ax.plot(x, y)
#The outer ellipse
x, y = np.meshgrid(np.arange(-1.0, 8.0, 0.025),
np.arange(-3.0, 6.5, 0.025))
ax.contour(x, y,
(Po[0][0] * x + Po[0][1] * y - co[0])**2 + (Po[1][0] * x + Po[1][1] * y - co[1])**2, [1])
ax.autoscale_view()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
fig.savefig('ellipsoid.png')
except:
print("Inner:")
print(" C = ", Ci)
print(" d = ", di)
print("Outer:")
print(" P = ", Po)
print(" c = ", co)
lpt.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: lpt.py
#
# Purpose: Demonstrates how to solve the multi-processor
# scheduling problem using the Fusion API.
##
import sys
import random
from mosek.fusion import *
def main():
#Parameters:
n = 30 #Number of tasks
m = 6 #Number of processors
lb = 1. #Range of short task lengths
ub = 5.
sh = 0.8 #Proportion of short tasks
n_short = int(n * sh)
n_long = n - n_short
random.seed(0)
T = sorted([random.uniform(lb, ub) for i in range(n_short)]
+ [random.uniform(20 * lb, 20 * ub) for i in range(n_long)], reverse=True)
print("# jobs(n) : ", n)
print("# machine(m): ", m)
with Model('Multi-processor scheduling') as M:
x = M.variable('x', [m, n], Domain.binary())
t = M.variable('t', 1, Domain.unbounded())
M.constraint(Expr.sum(x, 0), Domain.equalsTo(1.))
M.constraint(Expr.sub(Var.repeat(t, m), Expr.mul(x, T)),
Domain.greaterThan(0.))
M.objective(ObjectiveSense.Minimize, t)
#LPT heuristic
schedule = [0. for i in range(m)]
init = [0. for i in range(n * m)]
for i in range(n):
mm = schedule.index(min(schedule))
schedule[mm] += T[i]
init[n * mm + i] = 1.
#Comment this line to switch off feeding in the initial LPT solution
x.setLevel(init)
M.setLogHandler(sys.stdout)
M.setSolverParam("mioTolRelGap", .01)
M.solve()
print('initial solution:')
for i in range(m):
print('M', i, init[i * n:(i + 1) * n])
print('MOSEK solution:')
for i in range(m):
print('M', i, [y for y in x.slice([i, 0], [i + 1, n]).level()])
if __name__ == '__main__':
main()
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 sys
from mosek.fusion import *
import mosek.fusion.pythonic
with Model('mico1') as M:
x = M.variable(Domain.integral(Domain.unbounded()))
y = M.variable(Domain.integral(Domain.unbounded()))
t = M.variable()
M.constraint(Expr.vstack(t, x, y), Domain.inQCone())
M.constraint(Expr.vstack(x - 3.8, 1, y), Domain.inPExpCone())
M.objective(ObjectiveSense.Minimize, t)
M.setLogHandler(sys.stdout)
M.solve()
print('Solution: x = {0}, y = {1}'.format(x.level()[0], y.level()[0]))
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.
##
import sys
from mosek.fusion import *
def main(args):
A = [[50.0, 31.0],
[3.0, -2.0]]
c = [1.0, 0.64]
with Model('milo1') as M:
x = M.variable('x', 2, Domain.integral(Domain.greaterThan(0.0)))
# Create the constraints
# 50.0 x[0] + 31.0 x[1] <= 250.0
# 3.0 x[0] - 2.0 x[1] >= -4.0
M.constraint('c1', Expr.dot(A[0], x), Domain.lessThan(250.0))
M.constraint('c2', Expr.dot(A[1], x), Domain.greaterThan(-4.0))
# Set max solution time
M.setSolverParam('mioMaxTime', 60.0)
# Set max relative gap (to its default value)
M.setSolverParam('mioTolRelGap', 1e-4)
# Set max absolute gap (to its default value)
M.setSolverParam('mioTolAbsGap', 0.0)
# Set the objective function to (c^T * x)
M.objective('obj', ObjectiveSense.Maximize, Expr.dot(c, x))
# Solve the problem
M.solve()
# Get the solution values
print('[x0, x1] = ', x.level())
print("MIP rel gap = %.2f (%f)" % (M.getSolverDoubleInfo(
"mioObjRelGap"), M.getSolverDoubleInfo("mioObjAbsGap")))
if __name__ == '__main__':
main(sys.argv[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
# providing an initial feasible solution.
##
import sys
from mosek.fusion import *
def main(args):
c = [7.0, 10.0, 1.0, 5.0]
with Model('mioinitsol') as M:
n = 4
x = M.variable('x', n, Domain.greaterThan(0.0))
x.slice(0,3).makeInteger()
M.constraint(Expr.sum(x), Domain.lessThan(2.5))
# Set the objective function to (c^T * x)
M.objective('obj', ObjectiveSense.Maximize, Expr.dot(c, x))
# Assign values to integer variables.
# We only set a slice of xx
init_sol = [1.0, 1.0, 0.0]
x.slice(0,3).setLevel(init_sol)
# Request constructing the solution from integer variable values
M.setSolverParam("mioConstructSol", "on")
# Solve the problem
M.setLogHandler(sys.stdout)
M.solve()
# Get the solution values
ss = M.getPrimalSolutionStatus()
print("Solution status: {0}".format(ss))
sol = x.level()
print('x = {0}'.format(sol))
# Was the initial solution used?
constr = M.getSolverIntInfo("mioConstructSolution")
constrVal = M.getSolverDoubleInfo("mioConstructSolutionObj")
print("Construct solution utilization: {0}\nConstruct solution objective: {1:.3f}\n".format(constr, constrVal))
if __name__ == '__main__':
main(sys.argv[1:])
nearestcorr.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: nearestcorr.py
#
# Purpose:
# Solves the nearest correlation matrix problem
#
# minimize || A - X ||_F s.t. diag(X) = e, X is PSD
#
# as the equivalent conic program
#
# minimize t
#
# subject to (t, vec(A-X)) in Q
# diag(X) = e
# X >= 0
##
import sys
import mosek
import mosek.fusion
from mosek.fusion import *
import mosek.fusion.pythonic
from mosek import LinAlg
"""
Assuming that e is an NxN expression, return the lower triangular part as a vector.
"""
def vec(e):
N = e.getShape()[0]
msubi = range(N * (N + 1) // 2)
msubj = [i * N + j for i in range(N) for j in range(i + 1)]
mcof = [2.0**0.5 if i !=
j else 1.0 for i in range(N) for j in range(i + 1)]
S = Matrix.sparse(N * (N + 1) // 2, N * N, msubi, msubj, mcof)
return S @ Expr.flatten(e)
def nearestcorr(A):
N = A.numRows()
# Create a model
with Model("NearestCorrelation") as M:
# Setting up the variables
X = M.variable("X", Domain.inPSDCone(N))
t = M.variable("t", 1, Domain.unbounded())
# (t, vec (A-X)) \in Q
M.constraint("C1", Expr.vstack(t, vec(A-X)), Domain.inQCone())
# diag(X) = e
M.constraint("C2", X.diag() == 1.0)
# Objective: Minimize t
M.objective(ObjectiveSense.Minimize, t)
M.writeTask('nearcor.ptf')
M.solve()
return X.level(), t.level()
def nearestcorr_nucnorm(A, gammas):
N = A.numRows()
with Model("NucNorm") as M:
# Setup variables
t = M.variable("t", 1, Domain.unbounded())
X = M.variable("X", Domain.inPSDCone(N))
w = M.variable("w", N, Domain.greaterThan(0.0))
# D = diag(w)
D = Expr.mulElm(Matrix.eye(N), Var.repeat(w, 1, N))
# (t, vec (X + D - A)) in Q
M.constraint(Expr.vstack(t, vec(X+D-A)), Domain.inQCone())
result = []
for g in gammas:
# Objective: Minimize t + gamma*Tr(X)
M.objective(ObjectiveSense.Minimize, t + g * Expr.sum(X.diag()))
M.solve()
# Find eigenvalues of X and compute its rank
d = [0.0] * int(N)
LinAlg.syeig(mosek.uplo.lo, N, X.level(), d)
result.append(
(g, t.level()[0], sum([d[i] > 1e-6 for i in range(N)]), X.level()))
return result
if __name__ == '__main__':
N = 5
A = Matrix.dense(N, N, [0.0, 0.5, -0.1, -0.2, 0.5,
0.5, 1.25, -0.05, -0.1, 0.25,
-0.1, -0.05, 0.51, 0.02, -0.05,
-0.2, -0.1, 0.02, 0.54, -0.1,
0.5, 0.25, -0.05, -0.1, 1.25])
gammas = [0.1 * i for i in range(11)]
X, t = nearestcorr(A)
print("--- Nearest Correlation ---")
print("X = ")
print(X.reshape((N, N)))
print("t = ", t)
print("--- Nearest Correlation with Nuclear Norm---")
for g, res, rank, X in nearestcorr_nucnorm(A, gammas):
print("gamma=%f, res=%e, rank=%d" % (g, res, rank))
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
##
from mosek.fusion import *
import sys
if len(sys.argv) < 2:
print("Missing argument, syntax is:")
print(" python opt_server_sync.py addr [certpath]")
sys.exit(1)
serveraddr = sys.argv[1]
tlscert = None if len(sys.argv)<3 else sys.argv[2]
with Model('testOptServer') as M:
# Setup a simple test problem
x = M.variable('x', 3, Domain.greaterThan(0.0))
M.constraint("lc", Expr.dot([1.0, 1.0, 2.0], x), Domain.equalsTo(1.0))
M.objective("obj", ObjectiveSense.Minimize, Expr.sum(x))
# Attach log handler
M.setLogHandler(sys.stdout)
# Set OptServer URL
M.optserverHost(serveraddr)
# Path to certificate, if any
M.setSolverParam("remoteTlsCertPath", tlscert)
# Solve the problem on the OptServer
M.solve()
# Get the solution
print('x1,x2,x3 = %s' % str(x.level()))
parallel.py
#
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: parallel.py
#
# Purpose: Demonstrates parallel optimization using solveBatch()
#
from mosek.fusion import *
def makeToyParameterizedModel():
M = Model()
x = M.variable("x",3)
p = M.parameter("p")
M.objective(ObjectiveSense.Maximize, Expr.sum(x))
M.constraint(Expr.vstack(p,x), Domain.inQCone())
return M
# Example of how to use Model.solveBatch()
def main():
# Choose some sample parameters
n = 10 # Number of models to optimize
threadpoolsize = 4 # Total number of threads available
threadspermodel = 1 # Number of threads per each model
# Create a toy model for this example
M = makeToyParameterizedModel()
# Set up n copies of the model with different data
models = [M.clone() for _ in range(n)]
for i in range(n):
models[i].getParameter("p").setValue(i+1)
# We can set the number of threads individually per model
models[i].setSolverParam("numThreads", threadspermodel)
# Solve all models in parallel
status = Model.solveBatch(False, # No race
-1.0, # No time limit
threadpoolsize,
models) # Array of Models to solve
# Access information about each model
for i in range(n):
if status[i] == SolverStatus.OK:
print("Model {}: Status {}, Solution Status {}, Objective {:.3f}, Time {:.3f}".format(
i,
status[i],
models[i].getPrimalSolutionStatus(),
models[i].primalObjValue(),
models[i].getSolverDoubleInfo("optimizerTime")))
else:
print("Model {}: not solved".format(i))
if __name__ == "__main__":
main()
parameters.py
##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : parameters.py
#
# Purpose : Demonstrates a very simple example about how to set
# parameters and read information items
# with MOSEK Fusion
##
from mosek.fusion import *
# Create the Model
with Model() as M:
print('Test MOSEK parameter get/set functions')
# Set log level (integer parameter)
M.setSolverParam("log", 1)
# Select interior-point optimizer... (parameter with symbolic string values)
M.setSolverParam("optimizer", "intpnt")
# ... without basis identification (parameter with symbolic string values)
M.setSolverParam("intpntBasis", "never")
# Set relative gap tolerance (double parameter)
M.setSolverParam("intpntCoTolRelGap", 1.0e-7)
# The same in a different way
M.setSolverParam("intpntCoTolRelGap", "1.0e-7")
# Incorrect value
try:
M.setSolverParam("intpntCoTolRelGap", -1)
except ParameterError as e:
print('Wrong parameter value')
# Define and solve an optimization problem here
# M.solve()
# After optimization:
print('Get MOSEK information items')
tm = M.getSolverDoubleInfo("optimizerTime")
it = M.getSolverIntInfo("intpntIter")
print('Time: {0}\nIterations: {1}'.format(tm,it))
pinfeas.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: pinfeas.py
#
# Purpose: Demonstrates how to fetch a primal infeasibility certificate
# for a linear problem
#
##
from mosek.fusion import *
import sys
# Analyzes and prints infeasibility certificate for a single object,
# which can be a variable or constraint
def analyzeCertificate(name, # name of the analyzed object
size, # size of the object
duals, # double[], actual dual values
eps): # tolerance determining when a dual value is considered important
for i in range(0, size):
if abs(duals[i]) > eps:
print(f"{name}[{i}], dual = {duals[i]}")
# Construct the sample model from the example in the manual
sMat = Matrix.sparse(3, 7, [0,0,1,1,2,2,2],
[0,1,2,3,4,5,6],
[1.0]*7)
sBound = [200, 1000, 1000]
dMat = Matrix.sparse(4, 7, [0,0,1,2,2,3,3],
[0,4,1,2,5,3,6],
[1.0]*7)
dBound = [1100, 200, 500, 500]
c = [1, 2, 5, 2, 1, 2, 1]
M = Model("pinfeas")
x = M.variable("x", 7, Domain.greaterThan(0))
s = M.constraint("s", Expr.mul(sMat, x), Domain.lessThan(sBound))
d = M.constraint("d", Expr.mul(dMat, x), Domain.equalsTo(dBound))
M.objective(ObjectiveSense.Minimize, Expr.dot(c,x))
# Set useful debugging tools and solve the model
M.setLogHandler(sys.stdout)
M.writeTask("pinfeas.ptf")
M.solve()
# Check problem status
if M.getProblemStatus() == ProblemStatus.PrimalInfeasible:
# Set the tolerance at which we consider a dual value as essential
eps = 1e-7
# We want to retrieve infeasibility certificates
M.acceptedSolutionStatus(AccSolutionStatus.Certificate)
# Go through variable bounds
print("Variable bounds important for infeasibility: ")
analyzeCertificate(name = "x", size = x.getSize(), duals = x.dual(), eps = eps)
# Go through constraint bounds
print("Constraint bounds important for infeasibility: ")
analyzeCertificate(name = "s", size = s.getSize(), duals = s.dual(), eps = eps)
analyzeCertificate(name = "d", size = d.getSize(), duals = d.dual(), eps = eps)
else:
print("The problem is not primal infeasible, no certificate to show")
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.
#
##
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
"""
Purpose:
Computes the optimal portfolio for a given risk
Input:
n: Number of assets
mu: An n dimensional vector of expected returns
GT: A matrix with n columns so (GT')*GT = covariance matrix
x0: Initial holdings
w: Initial cash holding
gamma: Maximum risk (=std. dev) accepted
Output:
Optimal expected return and the optimal portfolio
"""
def BasicMarkowitz(n,mu,GT,x0,w,gamma):
with Model("Basic Markowitz") as M:
# Redirect log output from the solver to stdout for debugging.
# if uncommented.
# M.setLogHandler(sys.stdout)
# Defines the variables (holdings). Shortselling is not allowed.
x = M.variable("x", n, Domain.greaterThan(0.0))
# Maximize expected return
M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)
# The amount invested must be identical to initial wealth
M.constraint('budget', Expr.sum(x) == w+sum(x0))
# Imposes a bound on the risk
M.constraint('risk', Expr.vstack(gamma, GT @ x) == Domain.inQCone())
#M.constraint('risk', Expr.vstack(gamma, 0.5, Expr.mul(GT, x)), Domain.inRotatedQCone())
# Solves the model.
M.solve()
# Check if the solution is an optimal point
solsta = M.getPrimalSolutionStatus()
if (solsta != SolutionStatus.Optimal):
# See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
return np.dot(mu, x.level()), x.level()
if __name__ == '__main__':
n = 8
w = 59
mu = [0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628]
x0 = [8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0]
gammas = [36]
GT = [
[0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638],
[0. , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506],
[0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914],
[0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742],
[0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 ],
[0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187],
[0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327],
[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 ]
]
print("\n-----------------------------------------------------------------------------------")
print('Basic Markowitz portfolio optimization')
print("-----------------------------------------------------------------------------------\n")
for gamma in gammas:
er, x = BasicMarkowitz(n,mu,GT,x0,w,gamma)
print('Expected return: %.4e Std. deviation: %.4e' % (er, gamma))
print('Optimal portfolio: [%.4e, %.4e, %.4e, %.4e, %.4e, %.4e, %.4e, %.4e]' % tuple(x))
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.
#
##
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
import sys
"""
Purpose:
Computes several portfolios on the optimal portfolios by
for alpha in alphas:
maximize expected return - alpha * variance
subject to the constraints
Input:
n: Number of assets
mu: An n dimensional vector of expected returns
GT: A matrix with n columns so (GT')*GT = covariance matrix
x0: Initial holdings
w: Initial cash holding
alphas: List of the alphas
Output:
The efficient frontier as list of tuples (alpha, expected return, variance)
"""
def EfficientFrontier(n,mu,GT,x0,w,alphas):
with Model("Efficient frontier") as M:
frontier = []
# Defines the variables (holdings). Shortselling is not allowed.
x = M.variable("x", n, Domain.greaterThan(0.0)) # Portfolio variables
s = M.variable("s", 1, Domain.unbounded()) # Variance variable
# Total budget constraint
M.constraint('budget', Expr.sum(x) == w+sum(x0))
# Computes the risk
M.constraint('variance', Expr.vstack(s, 0.5, GT @ x), Domain.inRotatedQCone())
# Define objective as a weighted combination of return and variance
alpha = M.parameter()
M.objective('obj', ObjectiveSense.Maximize, x.T @ mu - s * alpha)
# Solve multiple instances by varying the parameter alpha
for a in alphas:
alpha.setValue(a)
M.solve()
# Check if the solution is an optimal point
solsta = M.getPrimalSolutionStatus()
if (solsta != SolutionStatus.Optimal):
# See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
frontier.append((a, np.dot(mu,x.level()), s.level()[0]))
return frontier
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 ]
]
# Some predefined alphas are chosen
alphas = [0.0, 0.01, 0.1, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 10.0]
frontier= EfficientFrontier(n,mu,GT,x0,w,alphas)
print("\n-----------------------------------------------------------------------------------")
print('Efficient frontier')
print("-----------------------------------------------------------------------------------\n")
print('%-12s %-12s %-12s' % ('alpha','return','risk (std. dev.)'))
for i in frontier:
print('%-12.4f %-12.4e %-12.4e' % (i[0],i[1],np.sqrt(i[2])))
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.
##
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
"""
Description:
Extends the basic Markowitz model with a market cost term.
Input:
n: Number of assets
mu: An n dimensional vector of expected returns
GT: A matrix with n columns so (GT')*GT = covariance matrix
x0: Initial holdings
w: Initial cash holding
gamma: Maximum risk (=std. dev) accepted
m: It is assumed that market impact cost for the j'th asset is
m_j|x_j-x0_j|^3/2
Output:
Optimal expected return and the optimal portfolio
"""
def MarkowitzWithMarketImpact(n,mu,GT,x0,w,gamma,m):
with Model("Markowitz portfolio with market impact") as M:
#M.setLogHandler(sys.stdout)
# Defines the variables. No shortselling is allowed.
x = M.variable("x", n, Domain.greaterThan(0.0))
# Variables computing market impact
t = M.variable("t", n, Domain.unbounded())
# Maximize expected return
M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)
# Invested amount + slippage cost = initial wealth
M.constraint('budget', Expr.sum(x) + Expr.dot(m,t) == w+sum(x0))
# Imposes a bound on the risk
M.constraint('risk', Expr.vstack(gamma, GT @ x), Domain.inQCone())
# t >= |x-x0|^1.5 using a power cone
M.constraint('tz', Expr.hstack(t, Expr.constTerm(n, 1.0), x-x0), Domain.inPPowerCone(2.0/3.0))
M.solve()
# Check if the solution is an optimal point
solsta = M.getPrimalSolutionStatus()
if (solsta != SolutionStatus.Optimal):
# See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
return x.level(), t.level()
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 ]
]
# Somewhat arbitrary choice of m
gamma = 0.36
m = n * [0.01]
xsol, tsol = MarkowitzWithMarketImpact(n,mu,GT,x0,w,gamma,m)
print("\n-----------------------------------------------------------------------------------")
print('Markowitz portfolio optimization with market impact cost')
print("-----------------------------------------------------------------------------------\n")
print('Expected return: %.4e Std. deviation: %.4e Market impact cost: %.4e' % \
(np.dot(mu,xsol),gamma,np.dot(m,tsol)))
print('Optimal portfolio: [%.4e, %.4e, %.4e, %.4e, %.4e, %.4e, %.4e, %.4e]' % tuple(xsol))
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.
##
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
"""
Description:
Extends the basic Markowitz model with a market cost term.
Input:
n: Number of assets
mu: An n dimensional vector of expected returns
GT: A matrix with n columns so (GT')*GT = covariance matrix
x0: Initial holdings
w: Initial cash holding
gamma: Maximum risk (=std. dev) accepted
f: If asset j is traded then a fixed cost f_j must be paid
g: If asset j is traded then a cost g_j must be paid for each unit traded
Output:
Optimal expected return and the optimal portfolio
"""
def MarkowitzWithTransactionsCost(n,mu,GT,x0,w,gamma,f,g):
# Upper bound on the traded amount
w0 = w+sum(x0)
u = n*[w0]
with Model("Markowitz portfolio with transaction costs") as M:
#M.setLogHandler(sys.stdout)
# Defines the variables. No shortselling is allowed.
x = M.variable("x", n, Domain.greaterThan(0.0))
# Additional "helper" variables
z = M.variable("z", n, Domain.unbounded())
# Binary variables
y = M.variable("y", n, Domain.binary())
# Maximize expected return
M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)
# Invest amount + transactions costs = initial wealth
M.constraint('budget', Expr.sum(x) + Expr.dot(f,y) + Expr.dot(g,z) == w0)
# Imposes a bound on the risk
M.constraint('risk', Expr.vstack(gamma, GT @ x), Domain.inQCone())
# z >= |x-x0|
M.constraint('buy', z >= x-x0)
M.constraint('sell', z >= x0-x)
# Alternatively, formulate the two constraints as
#M.constraint('trade', Expr.hstack(z, x-x0), Domain.inQcone())
# Constraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
M.constraint('y_on_off', z <= Expr.mulElm(u, y))
# Integer optimization problems can be very hard to solve so limiting the
# maximum amount of time is a valuable safe guard
M.setSolverParam('mioMaxTime', 180.0)
M.solve()
# Check if the solution is an optimal point
solsta = M.getPrimalSolutionStatus()
if (solsta != SolutionStatus.Optimal):
# See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
return x.level(), y.level(), z.level()
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 ]
]
f = n*[0.01]
g = n*[0.001]
gamma = 0.36
xsol, ysol, zsol = MarkowitzWithTransactionsCost(n,mu,GT,x0,w,gamma,f,g)
print("\n-----------------------------------------------------------------------------------")
print('Markowitz portfolio optimization with transactions cost')
print("-----------------------------------------------------------------------------------\n")
print('Expected return: %.4e Std. deviation: %.4e Transactions cost: %.4e' % \
(np.dot(mu,xsol),gamma, np.dot(f,ysol)+np.dot(g,zsol)))
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.
##
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
"""
Description:
Extends the basic Markowitz model with cardinality constraints.
Input:
n: Number of assets
mu: An n dimensional vector of expected returns
GT: A matrix with n columns so (GT')*GT = covariance matrix
x0: Initial holdings
w: Initial cash holding
gamma: Maximum risk (=std. dev) accepted
k: Maximum number of assets on which we allow to change position.
Output:
Optimal expected return and the optimal portfolio
"""
def MarkowitzWithCardinality(n,mu,GT,x0,w,gamma,KValues):
# Upper bound on the traded amount
w0 = w+sum(x0)
u = n*[w0]
with Model("Markowitz portfolio with cardinality bound") as M:
#M.setLogHandler(sys.stdout)
# Defines the variables. No shortselling is allowed.
x = M.variable("x", n, Domain.greaterThan(0.0))
# Additional "helper" variables
z = M.variable("z", n, Domain.unbounded())
# Binary variables - do we change position in assets
y = M.variable("y", n, Domain.binary())
# Maximize expected return
M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)
# The amount invested must be identical to initial wealth
M.constraint('budget', Expr.sum(x) == w+sum(x0))
# Imposes a bound on the risk
M.constraint('risk', Expr.vstack(gamma, GT @ x), Domain.inQCone())
# z >= |x-x0|
M.constraint('buy', z >= x-x0)
M.constraint('sell', z >= x0-x)
# Constraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
M.constraint('y_on_off', z <= Expr.mulElm(u,y))
# At most K assets change position
cardMax = M.parameter()
M.constraint('cardinality', Expr.sum(y)<= cardMax)
# Integer optimization problems can be very hard to solve so limiting the
# maximum amount of time is a valuable safe guard
M.setSolverParam('mioMaxTime', 180.0)
# Solve multiple instances by varying the parameter K
results = []
for K in KValues:
cardMax.setValue(K)
M.solve()
# Check if the solution is an optimal point
solsta = M.getPrimalSolutionStatus()
if (solsta != SolutionStatus.Optimal):
# See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
results.append(x.level())
return results
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.25
xsol = MarkowitzWithCardinality(n,mu,GT,x0,w,gamma,range(1,n+1))
print("\n-----------------------------------------------------------------------------------")
print('Markowitz portfolio optimization with cardinality constraints')
print("-----------------------------------------------------------------------------------\n")
for K in range(0,n):
print('Bound: {0} Expected return: {1:.4f} Solution {2}'.format(K+1, np.dot(mu, xsol[K]), xsol[K]))
portfolio_6_factor.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: portfolio_6_factor.py
#
# Description : Implements a basic portfolio optimization model
# with factor structured covariance matrix.
##
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
"""
Description:
Extends the basic Markowitz model with factor structure in the covariance matrix.
Input:
n: Number of securities
mu: An n dimensional vector of expected returns
G_factor_T: The factor (dense) part of the factorized risk
theta: specific risk vector
x0: Initial holdings
w: Initial cash holding
gamma: Maximum risk (=std. dev) accepted
Output:
Optimal expected return and the optimal portfolio
"""
def FactorModelMarkowitz(n, mu, G_factor_T, S_theta, x0, w, gamma):
with Model("Factor model Markowitz") as M:
# Variables
# The variable x is the fraction of holdings in each security.
# It is restricted to be positive, which imposes the constraint of no short-selling.
x = M.variable("x", n, Domain.greaterThan(0.0))
# Objective (quadratic utility version)
M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)
# Budget constraint
M.constraint('budget', Expr.sum(x) == w + sum(x0))
# Conic constraint for the portfolio std. dev
M.constraint('risk', Expr.vstack([Expr.constTerm(gamma),
G_factor_T @ x,
Expr.mulElm(np.sqrt(theta), x)]), Domain.inQCone())
# Solve optimization
M.solve()
# Check if the solution is an optimal point
solsta = M.getPrimalSolutionStatus()
if (solsta != SolutionStatus.Optimal):
# See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
return mu @ x.level(), x.level()
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
G_factor_T = G_factor.T
gammas = [0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48]
print("\n-----------------------------------------------------------------------------------")
print('Markowitz portfolio optimization with factor model')
print("-----------------------------------------------------------------------------------\n")
for gamma in gammas:
er, x = FactorModelMarkowitz(n, mu, G_factor_T, theta, x0, w, gamma)
print('Expected return: %.4e Std. deviation: %.4e' % (er, gamma))
np.set_printoptions(precision=4)
print(f'Optimal portfolio: {x}')
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
#
##
from mosek.fusion import *
with Model('pow1') as M:
x = M.variable('x', 3, Domain.unbounded())
x3 = M.variable()
x4 = M.variable()
# Create the linear constraint
M.constraint(Expr.dot(x, [1.0, 1.0, 0.5]), Domain.equalsTo(2.0))
# Create the power cone constraints
M.constraint(Var.vstack(x.slice(0,2), x3), Domain.inPPowerCone(0.2))
M.constraint(Expr.vstack(x.index(2), 1.0, x4), Domain.inPPowerCone(0.4))
# Set the objective function
M.objective(ObjectiveSense.Maximize, Expr.dot([1.0,1.0,-1.0], Var.vstack(x3, x4, x.index(0))))
# Solve the problem
M.solve()
# Get the linear solution values
solx = x.level()
print('x,y,z = %s' % str(solx))
primal_svm.py
#
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: primal_svm.py
#
# Purpose: Implements a simple soft-margin SVM
# using the Fusion API.
#
import mosek
from mosek.fusion import *
import mosek.fusion.pythonic
import random
def primal_svm(m,n,X,y,CC):
print("Number of data : %d"%m)
print("Number of features: %d"%n)
with Model() as M:
w = M.variable('w' , n, Domain.unbounded())
t = M.variable('t' , 1, Domain.unbounded())
b = M.variable('b' , 1, Domain.unbounded())
xi = M.variable('xi', m, Domain.greaterThan(0.))
M.constraint( Expr.mulElm( y, X @ w - Var.repeat(b,m) ) >= 1 - xi )
M.constraint( Expr.vstack(1., t, w), Domain.inRotatedQCone() )
print (' c | b | w ')
for C in CC:
M.objective( ObjectiveSense.Minimize, t + C * Expr.sum(xi) )
M.solve()
try:
cb = '{0:6} | {1:8f} | '.format(C,b.level()[0])
wstar =' '.join([ '{0:8f}'.format(wi) for wi in w.level()])
print (cb+wstar)
except:
pass;
if __name__ == '__main__':
CC=[ 500.0*i for i in range(10)]
m = 50
n = 3
seed= 0
random.seed(seed)
nump= random.randint(0,50)
numm= m - nump
y = [ 1. for i in range(nump)] + \
[ -1. for i in range(numm)]
mean = 1.
var = 1.
X= [ [ random.gauss( mean,var) for f in range(n) ] for i in range(nump)] + \
[ [ random.gauss(-mean,var) for f in range(n) ] for i in range(numm)]
primal_svm(m,n,X,y,CC)
qcqp_sdo_relaxation.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: qcqp_sdo_relaxation.py
#
# Purpose: Demonstrate how to use SDP to solve
# convex relaxation of a mixed-integer QCQP
##
import math
import sys
import numpy
import mosek
from mosek.fusion import *
import mosek.fusion.pythonic
# The relaxed SDP model
def miqcqp_sdo_relaxation(n,P,q):
M = Model()
M.setLogHandler(sys.stdout)
Z = M.variable("Z", Domain.inPSDCone(n+1))
X = Z[:-1, :-1]
x = Z[:-1, -1:]
M.constraint( X.diag() >= x )
M.constraint( Z[n, n] == 1 )
M.objective( ObjectiveSense.Minimize,
Expr.dot( P, X ) + 2.0 * Expr.dot(x, q) )
return M
# A direct integer model for minimizing |Ax-b|
def int_least_squares(n, A, b):
M = Model()
M.setLogHandler(sys.stdout)
x = M.variable("x", n, Domain.integral(Domain.unbounded()))
t = M.variable("t", 1, Domain.unbounded())
M.constraint( Expr.vstack(t, A @ x - b), Domain.inQCone() )
M.objective( ObjectiveSense.Minimize, t )
return M
# problem dimensions
n = 20
m = 2*n
# problem data
A = numpy.reshape(numpy.random.normal(0., 1.0, n*m), (m,n))
c = numpy.random.uniform(0., 1.0, n)
P = A.transpose().dot(A)
q = - P.dot(c)
b = A.dot(c)
# solve the problems
M = miqcqp_sdo_relaxation(n, P, q)
Mint = int_least_squares(n, A, b)
M.solve()
Mint.solve()
M.writeTask('M.ptf')
Mint.writeTask('Mint.ptf')
# rounded and optimal solution
xRound = numpy.rint(M.getVariable("Z").slice([0,n], [n,n+1]).level())
xOpt = numpy.rint(Mint.getVariable("x").level())
print(M.getSolverDoubleInfo("optimizerTime"), Mint.getSolverDoubleInfo("optimizerTime"))
print(numpy.linalg.norm(A.dot(xRound)-b), numpy.linalg.norm(A.dot(xOpt)-b))
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.
from mosek.fusion import *
# Problem data
c = [ 1.5, 2.5, 3.0 ]
A = [ [2, 4, 3],
[3, 2, 3],
[2, 3, 2] ]
b = [ 100000.0, 50000.0, 60000.0 ]
numvar = len(c)
numcon = len(b)
# Create a model and input data
with Model() as M:
x = M.variable("x",numvar, Domain.greaterThan(0.0))
con = M.constraint(Expr.mul(A, x), Domain.lessThan(b))
M.objective(ObjectiveSense.Maximize, Expr.dot(c, x))
# Solve the problem
M.solve()
print("x = {}".format(x.level()))
############# Replace an element of the A matrix ###############
x0 = x.index(0)
con.index(0).update(Expr.mul(3.0, x0), x0)
M.solve()
print("x = {}".format(x.level()))
############## Add a new variable ################
# Create a variable and a compound view of all variables
x3 = M.variable("y",Domain.greaterThan(0.0))
xNew = Var.vstack(x, x3)
print(x3.getShape())
# Add to the exising constraint
con.update(Expr.mul(x3, [4, 0, 1]), x3)
# Change the objective to include x3
M.objective(ObjectiveSense.Maximize, Expr.dot(c+[1.0], xNew))
M.solve()
print("x = {}".format(xNew.level()))
############## Add a new constraint ################
con2 = M.constraint(Expr.dot(xNew, [1, 2, 1, 1]), Domain.lessThan(30000.0))
M.solve()
print("x = {}".format(xNew.level()))
############## Change constraint bounds ################
cAll = Constraint.vstack(con, con2)
# Change bounds by effectively updating fixed terms with the difference
cAll.update([20000, 10000, 10000, 8000])
M.solve()
print("x = {}".format(xNew.level()))
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
from mosek.fusion import *
import sys
# Set up a small artificial conic problem to test with
def setupExample(M):
x = M.variable('x', 3, Domain.greaterThan(0.0))
y = M.variable('y', 3, Domain.unbounded())
z1 = Var.vstack(y.index(0), x.slice(0, 2))
z2 = Var.vstack(y.slice(1, 3), x.index(2))
M.constraint("lc", Expr.dot([1.0, 1.0, 2.0], x), Domain.equalsTo(1.0))
qc1 = M.constraint("qc1", z1, Domain.inQCone())
qc2 = M.constraint("qc2", z2, Domain.inRotatedQCone())
M.objective("obj", ObjectiveSense.Minimize, Expr.sum(y))
with Model() as M:
# (Optional) set a log stream
# M.setLoghandler(sys.stdout)
# (Optional) uncomment to see what happens when solution status is unknown
# M.setSolverParam("intpntMaxIterations", 1)
# In this example we set up a small conic problem
setupExample(M)
# Optimize
try:
M.solve()
# We expect solution status OPTIMAL (this is also default)
M.acceptedSolutionStatus(AccSolutionStatus.Optimal)
print("Optimal solution for x: {0}".format(M.getVariable('x').level()))
print("Optimal primal objective: {0}".format(M.primalObjValue()))
# .. continue analyzing the solution
except OptimizeError as e:
print("Optimization failed. Error: {0}".format(e))
except SolutionError as e:
# The solution with at least the expected status was not available.
# We try to diagnoze why.
print("Requested solution was not available.")
prosta = M.getProblemStatus()
if prosta == ProblemStatus.DualInfeasible:
print("Dual infeasibility certificate found.")
elif prosta == ProblemStatus.PrimalInfeasible:
print("Primal infeasibility certificate found.")
elif prosta == ProblemStatus.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(mosek.rescode(int(M.getSolverIntInfo("optimizeResponse"))))
print(" Termination code: {0} {1}".format(symname, desc))
else:
print("Another unexpected problem status {0} is obtained.".format(prosta))
except Exception as e:
print("Unexpected error: {0}".format(e))
sdo1.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: sdo1.py
#
# Purpose:
# Solves the mixed semidefinite and conic quadratic optimization problem
#
# [2, 1, 0]
# minimize Tr [1, 2, 1] * X + x0
# [0, 1, 2]
#
# [1, 0, 0]
# subject to Tr [0, 1, 0] * X + x0 = 1.0
# [0, 0, 1]
#
# [1, 1, 1]
# Tr [1, 1, 1] * X + x1 + x2 = 0.5
# [1, 1, 1]
#
# X >> 0, x0 >= (x1^2 + x2^2) ^ (1/2)
##
import sys
import mosek
from mosek.fusion import *
def main(args):
with Model("sdo1") as M:
# Setting up the variables
X = M.variable("X", Domain.inPSDCone(3))
x = M.variable("x", Domain.inQCone(3))
# Setting up constant coefficient matrices
C = Matrix.dense ( [[2.,1.,0.],[1.,2.,1.],[0.,1.,2.]] )
A1 = Matrix.eye(3)
A2 = Matrix.ones(3,3)
# Objective
M.objective(ObjectiveSense.Minimize, Expr.add(Expr.dot(C, X), x.index(0)))
# Constraints
M.constraint("c1", Expr.add(Expr.dot(A1, X), x.index(0)), Domain.equalsTo(1.0))
M.constraint("c2", Expr.add(Expr.dot(A2, X), Expr.sum(x.slice(1,3))), Domain.equalsTo(0.5))
M.solve()
print(X.level())
print(x.level())
if __name__ == '__main__':
main(sys.argv[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.fusion import *
import sys
import numpy as np
# Sample data in sparse, symmetric triplet format
C1_k = [0, 2]
C1_l = [0, 2]
C1_v = [1, 6]
A1_k = [0, 2, 0, 2]
A1_l = [0, 0, 2, 2]
A1_v = [1, 1, 1, 2]
C2_k = [0, 1, 0, 1, 2]
C2_l = [0, 0, 1, 1, 2]
C2_v = [1, -3, -3, 2, 1]
A2_k = [1, 0, 1, 3]
A2_l = [0, 1, 1, 3]
A2_v = [1, 1, -1, -3]
b = 23
k = -3
# Convert input data into Fusion sparse matrices
C1 = Matrix.sparse(3, 3, C1_k, C1_l, C1_v)
C2 = Matrix.sparse(4, 4, C2_k, C2_l, C2_v)
A1 = Matrix.sparse(3, 3, A1_k, A1_l, A1_v)
A2 = Matrix.sparse(4, 4, A2_k, A2_l, A2_v)
# Define the model
with Model('sdo2') as M:
# Two semidefinite variables
X1 = M.variable(Domain.inPSDCone(3))
X2 = M.variable(Domain.inPSDCone(4))
# Objective
M.objective(ObjectiveSense.Minimize, Expr.add(Expr.dot(C1,X1), Expr.dot(C2,X2)))
# Equality constraint
M.constraint(Expr.add([Expr.dot(A1,X1), Expr.dot(A2,X2)]), Domain.equalsTo(b))
# Inequality constraint
M.constraint(X2.index([0,1]), Domain.lessThan(k))
# Solve
M.setLogHandler(sys.stdout)
M.solve()
# Retrieve result
print("X1:\n{0}".format(np.reshape(X1.level(), (3,3))))
print("X2:\n{0}".format(np.reshape(X2.level(), (4,4))))
sdo3.py
#
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File : sdo3.py
#
# Purpose : Solves the semidefinite problem:
#
# min tr(X_1) + ... + tr(X_n)
# st. <A_11,X_1> + ... + <A_1n,X_n> >= b_1
# ...
# <A_k1,X_1> + ... + <A_kn,X_n> >= b_k
#
# where X_i are symmetric positive semidefinite of dimension d,
#
# A_ji are constant symmetric matrices and b_i are constant.
#
# This example is to demonstrate creating and using
# many matrix variables of the same dimension.
from mosek.fusion import *
import mosek.fusion.pythonic
import sys
import numpy as np
# Sample input data
n = 100
d = 4
k = 3
b = [9,10,11]
A = list(map(lambda x: x+np.transpose(x),
[np.random.normal(0, 5, size=(d, d)) for _ in range(k*n)]))
# Create a model with n semidefinite variables od dimension d x d
with Model("sdo3") as M:
X = M.variable(Domain.inPSDCone(d, n))
# Pick indexes of diagonal entries for the objective
alldiag = [[j,s,s] for j in range(n) for s in range(d)]
M.objective(ObjectiveSense.Minimize, Expr.sum( X.pick(alldiag) ))
# Each constraint is a sum of inner products
# Each semidefinite variable is a slice of X
for i in range(k):
M.constraint(Expr.add([Expr.dot(A[i*n+j], X[j:j+1, :, :].reshape([d,d]))
for j in range(n)]) >= b[i])
# Solve
M.setLogHandler(sys.stdout) # Add logging
M.writeTask("sdo3.ptf") # Save problem in readable format
M.solve()
# Get results. Each variable is a slice of X
print("Contributing variables:")
for j in range(n):
Xj = X.slice([j,0,0],[j+1,d,d]).level()
if any(Xj[s]>1e-6 for s in range(d*d)):
print("X{0}=\n{1}".format(j, np.reshape(Xj,(d,d))))
sospoly.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: sospoly.py
#
# Purpose:
# Models the cone of nonnegative polynomials and nonnegative trigonometric
# polynomials using Nesterov's framework [1].
#
# Given a set of coefficients (x0, x1, ..., xn), the functions model the
# cone of nonnegative polynomials
#
# P_m = { x \in R^{n+1} | x0 + x1*t + ... xn*t^n >= 0, forall t \in I }
#
# where I can be the entire real axis, the semi-infinite interval [0,inf), or
# a finite interval I = [a, b], respectively.
#
# References:
#
# [1] "Squared Functional Systems and Optimization Problems",
# Y. Nesterov, in High Performance Optimization,
# Kluwer Academic Publishers, 2000.
##
import mosek
from mosek.fusion import *
import mosek.fusion.pythonic
def Hankel(n, i, a=1.0):
'''Creates a Hankel matrix of dimension n+1, where
H_lk = a if l+k=i, and 0 otherwise.'''
if i < n+1:
return Matrix.sparse(n+1, n+1, range(i,-1,-1), range(i+1), (i+1)*[a])
return Matrix.sparse(n+1, n+1, range(n, i-n-1, -1), range(i-n, n+1), (2*n+1-i)*[a])
def nn_inf(M, x):
'''Models the cone of nonnegative polynomials on the real axis'''
m = int(x.getSize() - 1)
n = int(m/2) # degree of polynomial is 2*n
assert(m == 2*n)
# Setup variables
X = M.variable(Domain.inPSDCone(n+1))
# x_i = Tr H(n, i) * X i=0,...,m
for i in range(m+1):
M.constraint( x[i] == Expr.dot(Hankel(n,i), X) )
def nn_semiinf(M, x):
'''Models the cone of nonnegative polynomials on the semi-infinite interval [0,inf)'''
n = int(x.getSize()-1)
n1, n2 = int(n/2), int((n-1)/2)
# Setup variables
X1 = M.variable(Domain.inPSDCone(n1+1))
X2 = M.variable(Domain.inPSDCone(n2+1))
# x_i = Tr H(n1, i) * X1 + Tr H(n2,i-1) * X2, i=0,...,n
for i in range(n+1):
e1 = Expr.dot(Hankel(n1,i), X1)
e2 = Expr.dot(Hankel(n2,i-1),X2)
M.constraint( x[i] == e1 + e2 )
def nn_finite(M, x, a, b):
'''Models the cone of nonnegative polynomials on the finite interval [a,b]'''
assert(a < b)
m = int(x.getSize()-1)
n = int(m/2)
X1 = M.variable(Domain.inPSDCone(n+1))
if (m == 2*n):
X2 = M.variable(Domain.inPSDCone(n))
# x_i = Tr H(n,i)*X1 + (a+b)*Tr H(n-1,i-1) * X2 - a*b*Tr H(n-1,i)*X2 - Tr H(n-1,i-2)*X2, i=0,...,m
for i in range(m+1):
e1 = Expr.dot(Hankel(n , i ), X1) - Expr.dot(Hankel(n-1, i, a*b), X2)
e2 = Expr.dot(Hankel(n-1, i-1, a+b), X2) - Expr.dot(Hankel(n-1, i-2 ), X2)
M.constraint( x[i] == e1 + e2 )
else:
X2 = M.variable(Domain.inPSDCone(n+1))
# x_i = Tr H(n,i-1)*X1 - a*Tr H(n,i)*X1 + b*Tr H(n,i)*X2 - Tr H(n,i-1)*X2, i=0,...,m
for i in range(m+1):
e1 = Expr.dot(Hankel(n, i-1 ), X1) - Expr.dot(Hankel(n, i, a), X1)
e2 = Expr.dot(Hankel(n, i, b), X2) - Expr.dot(Hankel(n, i-1 ), X2)
M.constraint( x[i] == e1 + e2 )
def diff(M, x):
'''returns variables u representing the derivative of
x(0) + x(1)*t + ... + x(n)*t^n,
with u(0) = x(1), u(1) = 2*x(2), ..., u(n-1) = n*x(n).'''
n = int(x.getSize()-1)
u = M.variable(n, Domain.unbounded())
mask = Matrix.dense(n,1,[float(i) for i in range(1,n+1)])
M.constraint(u == Expr.mulElm(mask, x[1:n+1]))
return u
def fitpoly(data, n):
with Model('smooth poly') as M:
m = len(data)
A = [ [ data[j][0]**i for i in range(n+1) ] for j in range(m)]
b = [ data[j][1] for j in range(m) ]
x = M.variable('x', n+1, Domain.unbounded())
z = M.variable('z', 1, Domain.unbounded())
dx = diff(M, x)
M.constraint(A @ x == b)
# z - f'(t) >= 0, for all t \in [a, b]
ub = M.variable(n, Domain.unbounded())
dx0 = dx[0]
dx1n = dx[1:n]
M.constraint(ub == Expr.vstack(z - dx0, -dx1n))
nn_finite(M, ub, data[0][0], data[-1][0])
# f'(t) + z >= 0, for all t \in [a, b]
lb = M.variable(n, Domain.unbounded())
M.constraint(lb == Expr.vstack(z - dx0, dx1n))
nn_finite(M, lb, data[0][0], data[-1][0])
M.objective(ObjectiveSense.Minimize, z)
M.solve()
return x.level()
if __name__ == '__main__':
data = [ [-1.0, 1.0],
[ 0.0, 0.0],
[ 1.0, 1.0] ]
x2 = fitpoly(data, 2)
x4 = fitpoly(data, 4)
x8 = fitpoly(data, 8)
try:
from pyx import *
from pyx.graph.axis import painter, tick
text.set(mode="latex")
p = painter.regular(basepathattrs=[deco.earrow.normal],
titlepos=0.95,
titledist=-0.3,
titledirection=None,
outerticklength=graph.axis.painter.ticklength.normal,
)
g = graph.graphxy(width=8, xaxisat=0, yaxisat=0,
x=graph.axis.linear(title="$t$", min=-2.9, max=2.9,
painter=p),
y=graph.axis.linear(min=-0.05, max=1.9,
manualticks=[tick.tick(0, None, None),
tick.tick(0.5, label=r"$\frac{1}{2}$", ticklevel=0),
tick.tick(1.5, label=r"$\frac{3}{2}$", ticklevel=0),
],
painter=p))
def f(t,x): return t, sum([ x[i]*(t**i) for i in range(len(x)) ])
g.plot(graph.data.paramfunction("t", -3, 3, "x, y = f(t,x)", points=500, context={"f": f, "x":x2}),
[graph.style.line([style.linewidth.Thin, style.linestyle.solid])])
g.plot(graph.data.paramfunction("t", -3, 3, "x, y = f(t,x)", points=500, context={"f": f, "x":x4}),
[graph.style.line([style.linewidth.Thin, style.linestyle.solid])])
g.plot(graph.data.paramfunction("t", -3, 3, "x, y = f(t,x)", points=500, context={"f": f, "x":x8}),
[graph.style.line([style.linewidth.Thin, style.linestyle.solid])])
px, py = g.pos(1.3, f(1.3, x2)[1])
g.text(px - 0.1, py, "$f_2(t)$", [text.halign.right, text.valign.top])
px, py = g.pos(1.6, f(1.6, x4)[1])
g.text(px + 0.1, py, "$f_4(t)$", [text.halign.left, text.valign.top])
px, py = g.pos(1.31, f(1.31, x8)[1])
g.text(px - 0.1, py, "$f_8(t)$", [text.halign.right, text.valign.top])
g.writeEPSfile("sospoly")
g.writePDFfile("sospoly")
print("generated sospoly.eps")
except ImportError:
print("No PyX available")
# No pyx available
pass
sudoku.py
#
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: sudoku.py
#
# Purpose: A MILP-based SUDOKU solver
#
#
import sys
import mosek
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as npy
def print_solution(m, x):
n = m * m
print("\n")
for i in range(n):
ss = ""
for j in range(n):
if j % m == 0:
ss = ss + " |"
for k in range(n):
if x.index([i, j, k]).level() > 0.5:
ss = ss + " " + str(k + 1)
break
print(ss + ' |')
if (i + 1) % m == 0:
print("\n")
def main():
m = 3
n = m * m
hr_fixed = [[1, 5, 4],
[2, 2, 5], [2, 3, 8], [2, 6, 3],
[3, 2, 1], [3, 4, 2], [3, 5, 8], [3, 7, 9],
[4, 2, 7], [4, 3, 3], [4, 4, 1], [4, 7, 8], [4, 8, 4],
[6, 2, 4], [6, 3, 1], [6, 6, 9], [6, 7, 2], [6, 8, 7],
[7, 3, 4], [7, 5, 6], [7, 6, 5], [7, 8, 8],
[8, 4, 4], [8, 7, 1], [8, 8, 6],
[9, 5, 9]
]
fixed = [[f[0] - 1, f[1] - 1, f[2] - 1] for f in hr_fixed]
with Model('SUDOKU') as M:
x = M.variable("x",[n, n, n], Domain.binary())
#each value only once per dimension
for d in range(m):
M.constraint("row_sum(%d)" % d, Expr.sum(x, d) == 1)
#each number must appear only once in a block
for k in range(n):
for i in range(m):
for j in range(m):
M.constraint("blocksum(%d,%d,k=%d)" % (i,j,k),
Expr.sum(x[i*m : (i+1)*m, j*m : (j+1)*m, k : (k+1)]) == 1)
M.constraint("fix", x.pick(fixed) == 1.0)
M.setLogHandler(sys.stdout)
M.solve()
M.writeTask("sudoku.task")
#print the solution, if any...
if M.getPrimalSolutionStatus() in [SolutionStatus.Optimal]:
print_solution(m, x)
else:
print("No solution found!")
if __name__ == '__main__':
main()
total_variation.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: total_variation.py
#
# Purpose: Demonstrates how to solve a total
# variation problem using the Fusion API.
##
import sys
import mosek
from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np
def total_var(n,m):
M = Model('TV')
u = M.variable("u", [n+1,m+1], Domain.inRange(0.,1.0) )
t = M.variable("t", [n,m], Domain.unbounded() )
# In this example we define sigma and the input image f as parameters
# to demonstrate how to solve the same model with many data variants.
# Of course they could simply be passed as ordinary arrays if that is not needed.
sigma = M.parameter("sigma")
f = M.parameter("f", [n,m])
ucore = u[:-1, :-1]
deltax = u[1:, :-1] - ucore
deltay = u[:-1, 1:] - ucore
M.constraint( Expr.stack(2, t, deltax, deltay), Domain.inQCone().axis(2) )
M.constraint( Expr.vstack(sigma, Expr.flatten( f - ucore )),
Domain.inQCone() )
M.objective( ObjectiveSense.Minimize, Expr.sum(t) )
return M
#Display
def show(n,m,grid):
try:
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.imshow(grid, extent=(0,m,0,n),
interpolation='nearest', cmap=cm.jet,
vmin=0, vmax=1)
plt.show()
except:
print (grid)
if __name__ == '__main__':
np.random.seed(0)
n, m = 100, 200
# Create a parametrized model with given shape
M = total_var(n, m)
sigma = M.getParameter("sigma")
f = M.getParameter("f")
ucore = M.getVariable("u").slice([0,0], [n,m])
# Example: Linear signal with Gaussian noise
signal = np.reshape([[1.0*(i+j)/(n+m) for i in range(m)] for j in range(n)], (n,m))
noise = np.random.normal(0., 0.08, (n,m))
fVal = signal + noise
# Uncomment to get graphics:
# show(n, m, signal)
# show(n, m, fVal)
# Set value for f
f.setValue(fVal)
for sigmaVal in [0.0004, 0.0005, 0.0006]:
# Set new value for sigma and solve
sigma.setValue(sigmaVal*n*m)
M.solve()
sol = np.reshape(ucore.level(), (n,m))
# Now use the solution
# ...
# Uncomment to get graphics:
# show(n, m, np.reshape(ucore.level(), (n,m)))
print("rel_sigma = {sigmaVal} total_var = {var}".format(sigmaVal=sigmaVal,
var=M.primalObjValue() ))
M.dispose()
tsp.py
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File: tsp.py
#
# Purpose: Demonstrates a simple technique to the TSP
# usign the Fusion API.
#
##
from __future__ import print_function
import mosek
from mosek.fusion import *
import mosek.fusion.pythonic
import sys
def tsp(n, A, C, remove_selfloops, remove_2_hop_loops):
with Model() as M:
M.setLogHandler(sys.stdout)
x = M.variable([n,n], Domain.binary())
M.constraint(Expr.sum(x,0) == 1.0)
M.constraint(Expr.sum(x,1) == 1.0)
M.constraint(x <= A)
M.objective(ObjectiveSense.Minimize, Expr.dot(C, x))
if remove_2_hop_loops:
M.constraint(x + x.T <= 1.0)
if remove_selfloops:
M.constraint(x.diag() == 0)
it = 1
M.writeTask("tsp-0-%s-%s.ptf" % ('t' if remove_selfloops else 'f', 't' if remove_2_hop_loops else 'f'))
while True:
print("\n\n--------------------\nIteration",it)
M.solve()
print('\nsolution cost:', M.primalObjValue())
print('\nsolution:')
cycles = []
for i in range(n):
xi = x[i:i+1, :]
print(xi.level())
for j in range(n):
if xi.level()[j] <= 0.5 : continue
found = False
for c in cycles:
if len( [ a for a in c if i in a or j in a ] )> 0:
c.append( [i,j] )
found = True
break
if not found:
cycles.append([ [ i,j ]])
print('\ncycles:')
print([c for c in cycles])
if len(cycles)==1:
break;
for c in cycles:
M.constraint(Expr.sum(x.pick(c)) <= len(c) - 1)
it = it +1
return x.level(), c
return [],[]
def main():
A_i = [0,1,2,3,1,0,2,0]
A_j = [1,2,3,0,0,2,1,3]
C_v = [1.,1.,1.,1.,0.1,0.1,0.1,0.1]
n = max(max(A_i),max(A_j))+1
costs = Matrix.sparse(n,n,A_i,A_j,C_v)
x,c = tsp(n, Matrix.sparse(n,n,A_i,A_j,1.), costs , True, True)
x,c = tsp(n, Matrix.sparse(n,n,A_i,A_j,1.), costs , True, False)
if __name__ == '__main__':
main()
