16 List of examples

List of examples shipped in the distribution of Fusion API for Python:

Table 16.1 List of distributed examples
File Description
TrafficNetworkModel.py Demonstrates a traffic network problem as a conic quadratic problem (CQO)
alan.py A portfolio choice model alan.gms from the GAMS model library
baker.py A small bakery revenue maximization linear problem
breaksolver.py Shows how to break a long-running task
callback.py An example of data/progress callback
ceo1.py A simple conic exponential problem
cqo1.py A simple conic quadratic problem
diet.py Solving Stigler’s Nutrition model diet from the GAMS model library
duality.py Shows how to access the dual solution
elastic.py Linear regression with elastic net. Demonstrates model parametrization.
facility_location.py Demonstrates a small one-facility location problem (CQO)
gp1.py A simple geometric program (GP) in conic form
lo1.py A simple linear problem
logistic.py Implements logistic regression and simple log-sum-exp (CEO)
lownerjohn_ellipsoid.py Computes the Lowner-John inner and outer ellipsoidal approximations of a polytope (SDO, Power Cone)
lpt.py Demonstrates how to solve the multi-processor scheduling problem and input an integer feasible point (MIP)
mico1.py A simple mixed-integer conic problem
milo1.py A simple mixed-integer linear problem
mioinitsol.py A simple mixed-integer linear problem with an initial guess
modelLib.py Library of implementations of basic functions
nearestcorr.py Solves the nearest correlation matrix problem (SDO, CQO)
opt_server_sync.py Uses MOSEK OptServer to solve an optimization problem synchronously
parameters.py Shows how to set optimizer parameters and read information items
portfolio_1_basic.py Portfolio optimization - basic Markowitz model
portfolio_2_frontier.py Portfolio optimization - efficient frontier
portfolio_3_impact.py Portfolio optimization - market impact costs
portfolio_4_transcost.py Portfolio optimization - transaction costs
portfolio_5_card.py Portfolio optimization - cardinality constraints
pow1.py A simple power cone problem
primal_svm.py Implements a simple soft-margin Support Vector Machine (CQO)
qcqp_sdo_relaxation.py Demonstrate how to use SDP to solve convex relaxation of a mixed-integer QCQO problem
reoptimization.py Demonstrate how to modify and re-optimize a linear problem
response.py Demonstrates proper response handling
sdo1.py A simple semidefinite problem with one matrix variable and a quadratic cone
sdo2.py A simple semidefinite problem with two matrix variables
sdo3.py A simple semidefinite problem with many matrix variables of the same dimension
sospoly.py Models the cone of nonnegative polynomials and nonnegative trigonometric polynomials using Nesterov’s framework
sudoku.py A SUDOKU solver (MIP)
total_variation.py Demonstrates how to solve a total variation problem (CQO)
tsp.py 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

Listing 16.1 TrafficNetworkModel.py Click here to download.
##
# 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 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.mul(1.0 / T, Expr.add(Expr.dot(basetime, x), Expr.sum(d))))

        # 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)",
                    Expr.condense(
                       Expr.sub(
                         Expr.add(z, Expr.mulElm(x, cs_inv_matrix)),
                         s_inv_matrix)),
                     Domain.equalsTo(0.0))
        # Constraint (2)
        M.constraint("(2)",
                     Expr.sub( Expr.sum(x,0),
                               Expr.sum(x,1) ),
                     Domain.equalsTo(0.0))
        # Constraint (3)
        M.constraint("(3)", x.index(sink_idx, source_idx), Domain.equalsTo(T))

        M.writeTask("TrafficNetwork.opf")
        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

Listing 16.2 alan.py Click here to download.
##
# 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 *

#####################################################
### 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
M = Model('alan')
x = M.variable("x", numsec, Domain.greaterThan(0.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), Domain.equalsTo(1.0))
# define target expected return
mean = M.parameter(numsec)
M.constraint("dmean", Expr.dot(mean, x), Domain.greaterThan(target))

# Bound on risk
M.constraint("t > ||Ux||^2",
             Expr.vstack(0.5,
                         t,
                         Expr.mul(U, x)), Domain.inRotatedQCone())

##################################################
# 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)

# Finish
M.dispose()

baker.py

Listing 16.3 baker.py Click here to download.
# 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 *

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,
                    Expr.dot(revenue, production))

        # The production is constrained by stock:
        M.constraint(Expr.mul(recipe, production), Domain.lessThan(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

Listing 16.4 breaksolver.py Click here to download.
##
# 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

Listing 16.5 callback.py Click here to download.
##
#   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

Listing 16.6 ceo1.py Click here to download.
##
#   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 *

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.slice(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

Listing 16.7 cqo1.py Click here to download.
##
#   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 *

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.index(0), x.slice(0, 2))
    z2 = Var.vstack(y.slice(1, 3), x.index(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), Domain.equalsTo(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.opf')

    # 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

Listing 16.8 diet.py Click here to download.
# 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.
class DietModel(Model):

    def __init__(self, name,
                 foods,
                 nutrients,
                 daily_allowance,
                 nutritive_value):
        Model.__init__(self, name)
        finished = False
        try:
            self.foods = [str(f) for f in foods]
            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',
                                                 Set.make(self.foods),
                                                 Domain.greaterThan(0.0))
            self.__dailyNutrients = \
                self.constraint('Nutrient Balance',
                                Expr.mul(self.nutrientValue,
                                         self.__dailyPurchase),
                                Domain.greaterThan(self.dailyAllowance))
            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.]
    with DietModel('Stinglers Diet Problem',
                   foods,
                   nutrients,
                   daily_allowance,
                   nutritive_value) as M:
        M.solve()
        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))

duality.py

Listing 16.9 duality.py Click here to download.
# 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

Listing 16.10 elastic.py Click here to download.
# 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 *

# 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, Expr.sub(Expr.mul(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 = Expr.add([t, Expr.mul(lambda1, Expr.sum(p)), Expr.mul(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

Listing 16.11 facility_location.py Click here to download.
##
#  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

Listing 16.12 gp1.py Click here to download.
##
#   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 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), Domain.equalsTo(1.0))
    M.constraint(Expr.hstack(u,
                             Expr.constTerm(k, 1.0),
                             Expr.add(Expr.mul(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), Domain.lessThan(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

Listing 16.13 lo1.py Click here to download.
####
##  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 *

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.index(1), Domain.lessThan(10.0))
        M.constraint("c1", Expr.dot(A[0], x), Domain.equalsTo(30.0))
        M.constraint("c2", Expr.dot(A[1], x), Domain.greaterThan(15.0))
        M.constraint("c3", Expr.dot(A[2], x), Domain.lessThan(25.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
        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

Listing 16.14 logistic.py Click here to download.
# 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 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.constraint(Expr.add(z1, z2), Domain.equalsTo(1))
    M.constraint(Expr.hstack(z1, Expr.constTerm(n, 1.0), Expr.sub(u,t)), Domain.inPExpCone())
    M.constraint(Expr.hstack(z2, Expr.constTerm(n, 1.0), Expr.neg(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.add(Expr.sum(t), Expr.mul(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(Expr.mul(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

Listing 16.15 lownerjohn_ellipsoid.py Click here to download.
##
#  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 *

'''
Models the convex set

  S = { (x, t) \in R^n x R | x >= 0, t <= (x1 * x2 * ... * xn)^(1/n) }

using three-dimensional power cones
'''
def geometric_mean(M, x, t):
    n = int(x.getSize())
    if n==1:
      M.constraint(Expr.sub(t, x), Domain.lessThan(0.0))
    else:
      t2 = M.variable()
      M.constraint(Var.hstack(t2, x.index(n-1), t), Domain.inPPowerCone(1-1.0/n))
      geometric_mean(M, x.slice(0,n-1), t2)


'''
 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 number of rotated quadratic cones and affine hyperplanes,

   t <= (Z11*Z22*...*Znn)^{1/n}  (see geometric_mean).
'''
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.slice([0, 0], [n, n])
    Z   = Y.slice([0, n], [n, 2 * n])
    DZ  = Y.slice([n, n], [2 * 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(Expr.sub(DZ, Expr.mulElm(Z, Matrix.eye(n))), Domain.equalsTo(0.0))

    # t^n <= (Z11*Z22*...*Znn)
    geometric_mean(M, DZ.diag(), t)

    # 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(Expr.sub(b, Expr.mul(A, d)), Expr.mul(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),
                                 Expr.sub(Expr.mul(x, P),
                                          Expr.repeat(c,m,0))),
                     Domain.inQCone().axis(1))

        # 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 pentagon 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

Listing 16.16 lpt.py Click here to download.
##
#    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

Listing 16.17 mico1.py Click here to download.
##
#    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 *

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(Expr.sub(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

Listing 16.18 milo1.py Click here to download.
##
#    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

Listing 16.19 mioinitsol.py Click here to download.
##
#    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)

        # 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("Initial solution utilization: {0}\nInitial solution objective: {1:.3f}\n".format(constr, constrVal))                    

if __name__ == '__main__':
    main(sys.argv[1:])

modelLib.py

Listing 16.20 modelLib.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      modelLib.py
#
#  Purpose: Library of simple building blocks in Mosek Fusion.
##
from mosek.fusion import *

# Duplicate variables
# x = y
def dup(M, x, y):
    M.constraint(Expr.sub(x,y), Domain.equalsTo(0.0))

# Absolute value
# t >= |x|, where t, x have the same shape
def abs(M, t, x):
    M.constraint(Expr.add(t,x), Domain.greaterThan(0.0))
    M.constraint(Expr.sub(t,x), Domain.greaterThan(0.0))

# 1-norm
# t >= sum( |x_i| ), x is a vector expression
def norm1(M, t, x):
    u = M.variable(x.getShape(), Domain.unbounded())
    abs(M, u, x)
    M.constraint(Expr.sub(t, Expr.sum(u)), Domain.greaterThan(0.0))

# Square
# t >= x^2
def sq(M, t, x):
    M.constraint(Expr.hstack(0.5, t, x), Domain.inRotatedQCone())

# 2-norm
# t >= sqrt(x_1^2 + ... + x_n^2) where x is a vector
def norm2(M, t, x):
    M.constraint(Expr.vstack(t, x), Domain.inQCone())

# Power with exponent > 1
# t >= |x|^p (where p>1)
def pow(M, t, x, p):
    M.constraint(Expr.hstack(t, 1, x), Domain.inPPowerCone(1.0/p))

# Inverse of power 
# t >= 1/|x|^p, x>0 (where p>0)
def pow_inv(M, t, x, p):
    M.constraint(Expr.hstack(t, x, 1), Domain.inPPowerCone(1.0/(1.0+p)))

# p-norm, p>1
# t >= \|x\|_p (where p>1), x is a vector expression
def pnorm(M, t, x, p):
    n = int(x.getSize())
    r = M.variable(n)
    M.constraint(Expr.sub(t, Expr.sum(r)), Domain.equalsTo(0.0))
    M.constraint(Expr.hstack(Var.repeat(t,n), r, x), Domain.inPPowerCone(1.0-1.0/p))

# Geometric mean
# |t| <= (x_1...x_n)^(1/n), x_i>=0, x is a vector expression of length >= 1
def geo_mean(M, t, x):
    n = int(x.getSize())
    if n==1:
        abs(M, x, t)
    else:
        t2 = M.variable()
        M.constraint(Expr.hstack(t2, x.index(n-1), t), Domain.inPPowerCone(1.0-1.0/n))
        geo_mean(M, t2, x.slice(0,n-1))

# Logarithm
# t <= log(x), x>=0
def log(M, t, x):
    M.constraint(Expr.hstack(x, 1, t), Domain.inPExpCone())

# Exponential
# t >= exp(x)
def exp(M, t, x):
    M.constraint(Expr.hstack(t, 1, x), Domain.inPExpCone())

# Entropy
# t >= x * log(x), x>=0
def ent(M, t, x):
    M.constraint(Expr.hstack(1, x, Expr.neg(t)), Domain.inPExpCone())

# Relative entropy
# t >= x * log(x/y), x,y>=0
def relent(M, t, x, y):
    M.constraint(Expr.hstack(y, x, Expr.neg(t)), Domain.inPExpCone())

# Log-sum-exp
# log( sum_i(exp(x_i)) ) <= t, where x is a vector
def logsumexp(M, t, x):
    n = int(x.getSize())
    u = M.variable(n)
    M.constraint(Expr.hstack(u, Expr.constTerm(n, 1.0), Expr.sub(x, Var.repeat(t, n))), Domain.inPExpCone())
    M.constraint(Expr.sum(u), Domain.lessThan(1.0))

# Semicontinuous variable
# x = 0 or a <= x <= b
def semicontinuous(M, x, a, b):
    u = M.variable(x.getShape(), Domain.binary())
    M.constraint(Expr.sub(x, Expr.mul(a, u)), Domain.greaterThan(0.0))
    M.constraint(Expr.sub(x, Expr.mul(b, u)), Domain.lessThan(0.0))

# Indicator variable
# x!=0 implies t=1. Assumes that |x|<=1 in advance.
def indicator(M, t, x):
    M.constraint(t, Domain.inRange(0,1))
    t.makeInteger()
    abs(M, t, x)

# Logical OR
# x OR y, where x, y are binary
def logic_or(M, x, y):
    M.constraint(Expr.add(x, y), Domain.greaterThan(1.0))
# x_1 OR ... OR x_n, where x is a binary vector
def logic_or_vect(M, x):
    M.constraint(Expr.sum(x), Domain.greaterThan(1.0))

# SOS1 (NAND)
# at most one of x_1,...,x_n, where x is a binary vector (SOS1 constraint)
def logic_sos1(M, x):
    M.constraint(Expr.sum(x), Domain.lessThan(1.0))
# NOT(x AND y), where x, y are binary
def logic_nand(M, x, y):
    M.constraint(Expr.add(x, y), Domain.lessThan(1.0))

# Cardinality bound
# At most k of entries in x are nonzero, assuming in advance |x_i|<=1.
def card(M, x, k):
    t = M.variable(x.getShape(), Domain.binary())
    abs(M, t, x)
    M.constraint(Expr.sum(t), Domain.lessThan(k))

# This is just a syntactic test without much sense
def testModels():
    M = Model()
    x = M.variable()
    y = M.variable()
    t = M.variable()
    p = M.variable(5)
    a = M.variable([10,2])
    b = M.variable([10,2])
    e = M.variable(Domain.binary())
    f = M.variable(Domain.binary())

    log(M, t, x)
    exp(M, t, x)
    ent(M, t, x)
    relent(M, t, x, y)
    logsumexp(M, t, p)
    abs(M, a, b)
    norm1(M, t, a)
    sq(M, t, x)
    norm2(M, t, p)
    pow(M, t, x, 1.5)
    pow_inv(M, t, x, 3.3)
    geo_mean(M, t, p)
    semicontinuous(M, y, 1.1, 2.2)
    indicator(M, e, y)
    logic_or(M, e, f)
    logic_nand(M, e, f)
    card(M, b, 5)

# A sample problem using functions from the library
#
# max -sqrt(x^2 + y^2) + log(y) - x^1.5
#  st x >= y + 3
#
def testExample():
    M = Model()
    x = M.variable()
    y = M.variable()
    t = M.variable(3)

    M.constraint(Expr.sub(x, y), Domain.greaterThan(3.0))
    norm2(M, t.index(0), Var.vstack(x,y))
    log  (M, t.index(1), y)
    pow  (M, t.index(2), x, 1.5)

    M.objective(ObjectiveSense.Maximize, Expr.dot(t, [-1,1,-1]))

    import sys, numpy
    M.setLogHandler(sys.stdout)
    M.solve()
    print("x={0} y={1} obj={2}".format(x.level()[0], y.level()[0], M.primalObjValue()))
    if numpy.abs(M.primalObjValue()+10.5033852146)>1e-4:
        print("Wrong answer")
        sys.exit(1)
    else:
        print("OK")

if __name__ == "__main__":
    testModels()
    testExample()

nearestcorr.py

Listing 16.21 nearestcorr.py Click here to download.
##
#  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 *
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 Expr.mul(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
        v = vec(Expr.sub(A, X))
        M.constraint("C1", Expr.vstack(t, v), Domain.inQCone())

        # diag(X) = e
        M.constraint("C2", X.diag(), Domain.equalsTo(1.0))

        # Objective: Minimize t
        M.objective(ObjectiveSense.Minimize, t)
        M.writeTask('nearcor.task')
        M.writeTask('nearcor.cbf')
        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(Expr.sub(Expr.add(X, D), A))),
                     Domain.inQCone())

        result = []
        for g in gammas:
            # Objective: Minimize t + gamma*Tr(X)
            M.objective(ObjectiveSense.Minimize, Expr.add(
                t, Expr.mul(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(), 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

Listing 16.22 opt_server_sync.py Click here to download.
##
#  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) < 1:
    print("Missing argument, syntax is:")
    print("  python opt_server_sync.py serveraddr")
    sys.exit(1)

addr = sys.argv[1]

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(addr)

    # Solve the problem on the OptServer
    M.solve()

    # Get the solution
    print('x1,x2,x3 = %s' % str(x.level()))

parameters.py

Listing 16.23 parameters.py Click here to download.
##
#  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))

portfolio_1_basic.py

Listing 16.24 portfolio_1_basic.py Click here to download.
##
# 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 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, Expr.dot(mu,x))
        
        # The amount invested  must be identical to initial wealth
        M.constraint('budget', Expr.sum(x), Domain.equalsTo(w+sum(x0)))
        
        # Imposes a bound on the risk
        M.constraint('risk', Expr.vstack( gamma,Expr.mul(GT,x)), Domain.inQCone())

        # Solves the model.
        M.solve()

        return np.dot(mu,x.level())


if __name__ == '__main__':    

    n      = 3;
    w      = 1.0;   
    mu     = [0.1073,0.0737,0.0627]
    x0     = [0.0,0.0,0.0]
    gammas = [0.035,0.040,0.050,0.060,0.070,0.080,0.090]
    GT     = [
        [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 ],
        [ 0.0              , 0.102863378954911  , -0.00222873156550421],
        [ 0.0              , 0.0                ,  0.0338148677744977 ]
    ]


    print("\n-----------------------------------------------------------------------------------");
    print('Basic Markowitz portfolio optimization')
    print("-----------------------------------------------------------------------------------\n");
    for gamma in gammas:
        er = BasicMarkowitz(n,mu,GT,x0,w,gamma)
        print('Expected return: %.4e Std. deviation: %.4e' % (er,gamma))

portfolio_2_frontier.py

Listing 16.25 portfolio_2_frontier.py Click here to download.
##
# 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 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), Domain.equalsTo(w+sum(x0)))

        # Computes the risk
        M.constraint('variance', Expr.vstack(s, 0.5, Expr.mul(GT,x)), Domain.inRotatedQCone())

        # Define objective as a weighted combination of return and variance
        alpha = M.parameter()
        M.objective('obj', ObjectiveSense.Maximize, Expr.sub(Expr.dot(mu,x), Expr.mul(alpha,s)))
        
        # Solve multiple instances by varying the parameter alpha
        for a in alphas:
            alpha.setValue(a);

            M.solve()

            frontier.append((a, np.dot(mu,x.level()), s.level()[0]))

        return frontier


if __name__ == '__main__':

    n      = 3;
    w      = 1.0;
    mu     = [0.1073,0.0737,0.0627]
    x0     = [0.0,0.0,0.0]
    GT     = [
        [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 ],
        [ 0.0              , 0.102863378954911  , -0.00222873156550421],
        [ 0.0              , 0.0                ,  0.0338148677744977 ]
    ]

    # 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'))
    for i in frontier:
        print('%-12.4f  %-12.4e  %-12.4e' % (i[0],i[1],i[2]))

portfolio_3_impact.py

Listing 16.26 portfolio_3_impact.py Click here to download.
##
# 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 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, Expr.dot(mu,x))

        # Invested amount + slippage cost = initial wealth
        M.constraint('budget', Expr.add(Expr.sum(x),Expr.dot(m,t)), Domain.equalsTo(w+sum(x0)))

        # Imposes a bound on the risk
        M.constraint('risk', Expr.vstack(gamma,Expr.mul(GT,x)), Domain.inQCone())

        # t >= |x-x0|^1.5 using a power cone
        M.constraint('tz', Expr.hstack(t, Expr.constTerm(n, 1.0), Expr.sub(x,x0)), Domain.inPPowerCone(2.0/3.0))

        M.solve()

        return x.level(), t.level()


if __name__ == '__main__':    

    n      = 3;
    w      = 1.0;   
    mu     = [0.1073,0.0737,0.0627]
    x0     = [0.0,0.0,0.0]
    gammas = [0.035,0.040,0.050,0.060,0.070,0.080,0.090]
    GT     = [
        [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 ],
        [ 0.0              , 0.102863378954911  , -0.00222873156550421],
        [ 0.0              , 0.0                ,  0.0338148677744977 ]
    ]
                  
    # Somewhat arbitrary choice of m
    gamma = gammas[0]
    m = n*[1.0e-2]
    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)))

portfolio_4_transcost.py

Listing 16.27 portfolio_4_transcost.py Click here to download.
##
# 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 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, Expr.dot(mu,x))

        # Invest amount + transactions costs = initial wealth
        M.constraint('budget', Expr.add([ Expr.sum(x), Expr.dot(f,y),Expr.dot(g,z)] ), Domain.equalsTo(w0))

        # Imposes a bound on the risk
        M.constraint('risk', Expr.vstack( gamma,Expr.mul(GT,x)), Domain.inQCone())

        # z >= |x-x0| 
        M.constraint('buy', Expr.sub(z,Expr.sub(x,x0)),Domain.greaterThan(0.0))
        M.constraint('sell', Expr.sub(z,Expr.sub(x0,x)),Domain.greaterThan(0.0))
        # Alternatively, formulate the two constraints as
        #M.constraint('trade', Expr.hstack(z,Expr.sub(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', Expr.sub(z,Expr.mulElm(u,y)), Domain.lessThan(0.0))

        # 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()

        return x.level(), y.level(), z.level() 

if __name__ == '__main__':    

    n      = 3;
    w      = 1.0;   
    mu     = [0.1073,0.0737,0.0627]
    x0     = [0.0,0.0,0.0]
    gammas = [0.035,0.040,0.050,0.060,0.070,0.080,0.090]
    GT     = [
        [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 ],
        [ 0.0              , 0.102863378954911  , -0.00222873156550421],
        [ 0.0              , 0.0                ,  0.0338148677744977 ]
    ]


    f = n*[0.01]
    g = n*[0.001]
    gamma = gammas[0]
    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

Listing 16.28 portfolio_5_card.py Click here to download.
##
# 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 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, Expr.dot(mu,x))

        # The amount invested  must be identical to initial wealth
        M.constraint('budget', Expr.sum(x), Domain.equalsTo(w+sum(x0)))

        # Imposes a bound on the risk
        M.constraint('risk', Expr.vstack( gamma,Expr.mul(GT,x)), Domain.inQCone())

        # z >= |x-x0| 
        M.constraint('buy', Expr.sub(z,Expr.sub(x,x0)),Domain.greaterThan(0.0))
        M.constraint('sell', Expr.sub(z,Expr.sub(x0,x)),Domain.greaterThan(0.0))

        # 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', Expr.sub(z,Expr.mulElm(u,y)), Domain.lessThan(0.0))

        # At most k assets change position
        cardMax = M.parameter()
        M.constraint('cardinality', Expr.sub(Expr.sum(y), cardMax), Domain.lessThan(0))

        # 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()
            results.append(x.level())

        return results


if __name__ == '__main__':    

    n      = 3;
    w      = 1.0;   
    mu     = [0.1073,0.0737,0.0627]
    x0     = [0.0,0.0,0.0]
    gamma  = 0.035
    GT     = [
        [ 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 ],
        [ 0.0              , 0.102863378954911  , -0.00222873156550421],
        [ 0.0              , 0.0                ,  0.0338148677744977 ]
    ]

    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]))

pow1.py

Listing 16.29 pow1.py Click here to download.
##
#   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

Listing 16.30 primal_svm.py Click here to download.
#
# 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 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.add( 
                Expr.mulElm( y, 
                             Expr.sub( Expr.mul(X,w), Var.repeat(b,m) ) 
                         ),
                xi
            ), 
            Domain.greaterThan( 1. ) )

        M.constraint( Expr.vstack(1., t, w), Domain.inRotatedQCone() )

        print ('   c   |    b      |     w ')

        for C in CC:
            M.objective( ObjectiveSense.Minimize, Expr.add( t, Expr.mul(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

Listing 16.31 qcqp_sdo_relaxation.py Click here to download.
##
# 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 *

# 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.slice([0,0], [n,n])
    x = Z.slice([0,n], [n,n+1])

    M.constraint( Expr.sub(X.diag(), x), Domain.greaterThan(0.) )
    M.constraint( Z.index(n,n), Domain.equalsTo(1.) )

    M.objective( ObjectiveSense.Minimize, Expr.add( 
        Expr.sum( Expr.mulElm( P, X ) ), 
        Expr.mul( 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, Expr.sub(Expr.mul(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')
M.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

Listing 16.32 reoptimization.py Click here to download.
##
#  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)
        # 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

Listing 16.33 response.py Click here to download.
##
#  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

Listing 16.34 sdo1.py Click here to download.
##
#  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

Listing 16.35 sdo2.py Click here to download.
#
#  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

Listing 16.36 sdo3.py Click here to download.
#
#  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 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.slice([j,0,0], [j+1,d,d]).reshape([d,d])) 
                               for j in range(n)]), 
                     Domain.greaterThan(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

Listing 16.37 sospoly.py Click here to download.
##
#  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 *

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( Expr.sub(x.index(i), Expr.dot(Hankel(n,i),X)), Domain.equalsTo(0.0))

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( Expr.sub(x.index(i), Expr.add(e1, e2)), Domain.equalsTo(0.0) )
                 
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.sub(Expr.dot(Hankel(n  , i       ), X1), Expr.dot(Hankel(n-1, i, a*b), X2))
            e2 = Expr.sub(Expr.dot(Hankel(n-1, i-1, a+b), X2), Expr.dot(Hankel(n-1, i-2   ), X2))
            M.constraint( Expr.sub(x.index(i), Expr.add(e1, e2)), Domain.equalsTo(0.0) )
    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.sub(Expr.dot(Hankel(n, i-1 ), X1), Expr.dot(Hankel(n, i, a), X1))
            e2 = Expr.sub(Expr.dot(Hankel(n, i, b), X2), Expr.dot(Hankel(n, i-1 ), X2))
            M.constraint( Expr.sub(x.index(i), Expr.add(e1, e2)), Domain.equalsTo(0.0) )

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(Expr.sub(u, Expr.mulElm(mask, x.slice(1,n+1))), 
                 Domain.equalsTo(0.0))
    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(Expr.mul( A, x), Domain.equalsTo(b))
                        
        # z - f'(t) >= 0, for all t \in [a, b]
        ub = M.variable(n, Domain.unbounded())

        dx0= dx.index(0)
        dx1n= dx.slice(1,n)

        M.constraint(Expr.sub(ub, Expr.vstack(Expr.sub(z,dx0), Expr.neg(dx1n) ) ), 
                     Domain.equalsTo(0.0))

        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(Expr.sub(lb, Expr.vstack(Expr.add(z, dx0), dx1n)), 
                     Domain.equalsTo(0.0))

        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

Listing 16.38 sudoku.py Click here to download.
#
# 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 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), Domain.equalsTo(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.slice([i * m, j * m, k], [(i + 1) * m, (j + 1) * m, k + 1])),
                                 Domain.equalsTo(1.))

        M.constraint("fix",x.pick(fixed), Domain.equalsTo(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

Listing 16.39 total_variation.py Click here to download.
##
# 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 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.slice( [0,0], [n,m] )

    deltax = Expr.sub( u.slice( [1,0], [n+1,m] ), ucore)
    deltay = Expr.sub( u.slice( [0,1], [n,m+1] ), ucore)

    M.constraint( Expr.stack(2, t, deltax, deltay), Domain.inQCone().axis(2) )

    M.constraint( Expr.vstack(sigma, Expr.flatten( Expr.sub( 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

Listing 16.40 tsp.py Click here to download.
##
#  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 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), Domain.equalsTo(1.0))
        M.constraint(Expr.sum(x,1), Domain.equalsTo(1.0))
        M.constraint(x, Domain.lessThan( A ))

        M.objective(ObjectiveSense.Minimize, Expr.dot(C ,x))

        if remove_2_hop_loops:
            M.constraint(Expr.add(x, x.transpose()), Domain.lessThan(1.0))

        if remove_selfloops:
            M.constraint(x.diag(), Domain.equalsTo(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.slice([i,0],[i+1,n])
                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)), Domain.lessThan( 1.0 * 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()