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

djc1.py

A simple problem with disjunctive constraints (DJC)

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

nearestcorr.py

Solves the nearest correlation matrix problem (SDO, CQO)

opt_server_sync.py

Uses MOSEK OptServer to solve an optimization problem synchronously

parallel.py

Demonstrates parallel optimization using a batch method in MOSEK

parameters.py

Shows how to set optimizer parameters and read information items

pinfeas.py

Shows how to obtain and analyze a primal infeasibility certificate

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

portfolio_6_factor.py

Portfolio optimization - factor model

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 mosek.fusion.pythonic
import numpy as np
import sys


def main(args):
    n = 4
    arc_i = [0,    0,    2,    1,    2]
    arc_j = [1,    2,    1,    3,    3]
    arc_base = [4.0,  1.0,  2.0,  1.0,  6.0]
    arc_cap = [10.0, 12.0, 20.0, 15.0, 10.0]
    arc_sens = [0.1,  0.7,  0.9,  0.5,  0.1]

    T = 20.0
    source_idx = 0
    sink_idx = 3

    with Model() as M:
        narcs = len(arc_i)

        NxN = Set.make(n, n)
        sens = Matrix.sparse(n, n, arc_i, arc_j, arc_sens)
        cap = Matrix.sparse(n, n, arc_i, arc_j, arc_cap)
        basetime = Matrix.sparse(n, n, arc_i, arc_j, arc_base)

        cs_inv_matrix = \
            Matrix.sparse(n, n, arc_i, arc_j,
                          [1.0 / (arc_sens[i] * arc_cap[i]) for i in range(narcs)])
        s_inv_matrix = \
            Matrix.sparse(n, n, arc_i, arc_j,
                          [1.0 / arc_sens[i] for i in range(narcs)])

        sparsity = [ [ i, j] for i,j in zip(arc_i,arc_j) ]

        x = M.variable("traffic_flow", [n,n], Domain.greaterThan(0.0).sparse(sparsity+[ [sink_idx, source_idx] ]))
        t = M.variable("travel_time",  [n,n], Domain.greaterThan(0.0).sparse(sparsity))
        d = M.variable("d",            [n,n], Domain.greaterThan(0.0).sparse(sparsity))
        z = M.variable("z",            [n,n], Domain.greaterThan(0.0).sparse(sparsity))

        # Set the objective:
        M.objective("Average travel time",
                    ObjectiveSense.Minimize,
                    (Expr.dot(basetime, x) + Expr.sum(d)) / T)

        # Set up constraints
        # Constraint (1a)
        numnz = len(arc_sens)

        v = Var.hstack([ d.pick(sparsity),
                         z.pick(sparsity),
                         x.pick(sparsity) ])

        M.constraint("(1a)", v, Domain.inRotatedQCone(narcs,3))

        # Constraint (1b)
        M.constraint("(1b)", z + Expr.mulElm(x, cs_inv_matrix) ==  s_inv_matrix)

        # Constraint (2)
        M.constraint("(2)", Expr.sum(x,0) == Expr.sum(x,1))

        # Constraint (3)
        M.constraint("(3)", x[sink_idx, source_idx] == T)

        M.writeTask("TrafficNetwork.ptf")
        M.setLogHandler(sys.stdout);
        M.solve()
        M.writeTask("tn.ptf")

        M.writeTask("tn.ptf")
        flow = x.pick(sparsity).level()

        print("Optimal flow:")
        for i, j, v in zip(arc_i, arc_j, flow):
            print("\tflow node%d->node%d = %f" % (i, j, v))


main(sys.argv[1:])

alan.py

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 *
import mosek.fusion.pythonic

#####################################################
### Problem data:
# Security names
securities = ["hardware", "software", "show-biz", "t-bills"]
# Two examples of mean returns on securities
mean1 = [9.0, 7.0, 11.0, 5.0]
mean2 = [8.0, 9.0, 12.0, 7.0]
# Target mean return
target = 10.0
# Factor of covariance matrix.
U = Matrix.dense([[2.0, 1.5, -0.5, 0.0],
                  [0.0, 1.93649167, 0.90369611, 0.0],
                  [0.0, 0.0, 2.98886824, 0.0],
                  [0.0, 0.0, 0.0, 0.0]])

numsec = len(securities)

###################################################
# Create a model with expected returns as parameter
with Model('alan') as M:
    x = M.variable("x", numsec, Domain.greaterThan(0.0).withNamesOnAxis(securities,0))
    t = M.variable("t", 1, Domain.greaterThan(0.0))

    M.objective("minvar", ObjectiveSense.Minimize, t)

    # sum securities to 1.0
    M.constraint("wealth", Expr.sum(x) == 1.0)
    # define target expected return
    mean = M.parameter(numsec)
    M.constraint("dmean", Expr.dot(mean, x) >= target)

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

    M.writeTask('alan.ptf')
    ##################################################
    # Define a method that solves a single instance
    def solve(meanVal):
        print("Solve using mean={0}".format(meanVal))
        mean.setValue(meanVal)
        M.solve()
        solx = x.level()
        print("Solution: {0}".format(solx))
        return solx

    #################################################
    # Solve two instances with different means
    sol1 = solve(mean1)
    sol2 = solve(mean2)

baker.py

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 *
import mosek.fusion.pythonic

ingredientnames = ["Flour", "Sugar", "Butter"]
stock = [150.0, 22.0, 25.0]
recipe = [[3.0, 5.0],
          [1.0, 0.5],
          [1.2, 0.5]]
productnames = ["Cakes", "Breads"]
revenue = [4.0, 6.0]

def main(args):
    with Model("Recipe") as M:
        # "production" defines the amount of each product to bake.
        production = M.variable("production",
                                Set.make(productnames),
                                Domain.greaterThan(0.0))
        # The objective is to maximize the total revenue.
        M.objective("revenue",
                    ObjectiveSense.Maximize,
                    production.T @ revenue)

        # The production is constrained by stock:
        M.constraint(recipe @ production <= stock)
        # We set stdout as the loghandler
        M.setLogHandler(sys.stdout)
    # We solve and fetch the solution:
        M.solve()

        print("Solution:")
        res = production.level()
        for name, val in zip(productnames, res):
            print(" Number of %s : %.1f" % (name, val))
        print(" Revenue : %.2f" % (res[0] * revenue[0] + res[1] * revenue[1]))

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

breaksolver.py

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 *
import mosek.fusion.pythonic

with Model('ceo1') as M:

    x = M.variable('x', 3, Domain.unbounded())

    # Create the constraint
    #      x[0] + x[1] + x[2] = 1.0
    M.constraint("lc", Expr.sum(x), Domain.equalsTo(1.0))

    # Create the conic exponential constraint
    expc = M.constraint("expc", x, Domain.inPExpCone())

    # Set the objective function to (x[0] + x[1])
    M.objective("obj", ObjectiveSense.Minimize, Expr.sum(x[0:2]))

    # Solve the problem
    M.solve()

    M.writeTask("ceo1.ptf")

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

    # Get conic solution of expc
    expcval  = expc.level()
    expcdual = expc.dual()
    print('expc levels                = %s' % str(expcval))
    print('expc dual conic var levels = %s' % str(expcdual))

cqo1.py

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 *
import mosek.fusion.pythonic

with Model('cqo1') as M:

    x = M.variable('x', 3, Domain.greaterThan(0.0))
    y = M.variable('y', 3, Domain.unbounded())

    # Create the aliases
    #      z1 = [ y[0],x[0],x[1] ]
    #  and z2 = [ y[1],y[2],x[2] ]
    z1 = Var.vstack(y[0], x[0:2])
    z2 = Var.vstack(y[1:3], x[2])

    # Create the constraint
    #      x[0] + x[1] + 2.0 x[2] = 1.0
    M.constraint("lc", Expr.dot([1.0, 1.0, 2.0], x) == 1.0)

    # Create the constraints
    #      z1 belongs to C_3
    #      z2 belongs to K_3
    # where C_3 and K_3 are respectively the quadratic and
    # rotated quadratic cone of size 3, i.e.
    #                 z1[0] >= sqrt(z1[1]^2 + z1[2]^2)
    #  and  2.0 z2[0] z2[1] >= z2[2]^2
    qc1 = M.constraint("qc1", z1, Domain.inQCone())
    qc2 = M.constraint("qc2", z2, Domain.inRotatedQCone())

    # Set the objective function to (y[0] + y[1] + y[2])
    M.objective("obj", ObjectiveSense.Minimize, Expr.sum(y))

    # Solve the problem
    M.solve()
    M.writeTask('cqo1.ptf')

    # Get the linear solution values
    solx = x.level()
    soly = y.level()
    print('x1,x2,x3 = %s' % str(solx))
    print('y1,y2,y3 = %s' % str(soly))

    # Get conic solution of qc1
    qc1lvl = qc1.level()
    qc1sn = qc1.dual()
    print('qc1 levels                = %s' % str(qc1lvl))
    print('qc1 dual conic var levels = %s' % str(qc1sn))

diet.py

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.
def DietModel(name,
              foods,
              nutrients,
              daily_allowance,
              nutritive_value):
    with Model() as DM:
        dailyAllowance = numpy.array(daily_allowance, float)
        nutrientValue = Matrix.dense(nutritive_value).transpose()

        M = len(foods)
        N = len(nutrients)
        if len(dailyAllowance) != N:
            raise ValueError("Length of daily_allowance does not match "
                             "the number of nutrients")
        if nutrientValue.numColumns() != M:
            raise ValueError("Number of rows in nutrient_value does not "
                             "match the number of foods")
        if nutrientValue.numRows() != N:
            raise ValueError("Number of columns in nutrient_value does "
                             "not match the number of nutrients")

        dailyPurchase = DM.variable(
            'Daily Purchase',
            [M],
            Domain.greaterThan(0.0).withNamesOnAxis(foods,0))
        dailyNutrients = DM.constraint(
            'Nutrient Balance',
            Expr.mul(nutrientValue,
                     dailyPurchase),
            Domain.greaterThan(dailyAllowance).withNamesOnAxis(nutrients,0))
        DM.objective(ObjectiveSense.Minimize,
                    Expr.sum(dailyPurchase))

        DM.solve()
        DM.writeTask("diet.ptf")

        return dailyPurchase.level(),dailyNutrients.level()

# class DietModel(Model):

#     def __init__(self, name,
#                  foods,
#                  nutrients,
#                  daily_allowance,
#                  nutritive_value):
#         Model.__init__(self, name)
#         finished = False
#         try:
#             self.foods = numpy.array([str(f) for f in foods],dtype=numpy.dtype(object))
#             print(self.foods,len(numpy.array(self.foods,dtype=numpy.dtype(object))))
#             self.nutrients = [str(n) for n in nutrients]
#             self.dailyAllowance = numpy.array(daily_allowance, float)
#             self.nutrientValue = Matrix.dense(nutritive_value).transpose()

#             M = len(self.foods)
#             N = len(self.nutrients)
#             if len(self.dailyAllowance) != N:
#                 raise ValueError("Length of daily_allowance does not match "
#                                  "the number of nutrients")
#             if self.nutrientValue.numColumns() != M:
#                 raise ValueError("Number of rows in nutrient_value does not "
#                                  "match the number of foods")
#             if self.nutrientValue.numRows() != N:
#                 raise ValueError("Number of columns in nutrient_value does "
#                                  "not match the number of nutrients")

#             self.__dailyPurchase = self.variable('Daily Purchase',
#                                                  [M],
#                                                  Domain.greaterThan(0.0).withNamesOnAxis(self.foods,0))
#             self.__dailyNutrients = \
#                 self.constraint('Nutrient Balance',
#                                 Expr.mul(self.nutrientValue,
#                                          self.__dailyPurchase),
#                                 Domain.greaterThan(self.dailyAllowance).withNamesOnAxis(self.nutrients,0))
#             self.objective(ObjectiveSense.Minimize,
#                            Expr.sum(self.__dailyPurchase))
#             finished = True
#         finally:
#             if not finished:
#                 self.__del__()

#     def dailyPurchase(self):
#         return self.__dailyPurchase.level()

#     def dailyNutrients(self):
#         return self.__dailyNutrients.level()

if __name__ == '__main__':
    nutrient_unit = [
        "thousands", "grams", "grams",
        "milligrams", "thousand ius", "milligrams",
        "milligrams", "milligrams", "milligrams"]
    nutrients = [
        "calorie", "protein", "calcium",
        "iron", "vitamin a", "vitamin b1",
        "vitamin b2", "niacin", "vitamin c"]
    foods = [
        "wheat", "cornmeal", "cannedmilk", "margarine", "cheese",
        "peanut butter", "lard", "liver", "porkroast", "salmon",
        "greenbeans", "cabbage", "onions", "potatoes", "spinach",
        "sweet potatos", "peaches", "prunes", "limabeans", "navybeans"]
    nutritive_value = [
        #  calorie       calcium      vitamin a        vitamin b2      vitamin c
        #         protein        iron           vitamin b1      niacin
        [44.7, 1411, 2.0, 365, 0, 55.4, 33.3, 441, 0],  # wheat
        [36, 897, 1.7, 99, 30.9, 17.4, 7.9, 106, 0],  # cornmeal
        [8.4, 422, 15.1, 9, 26, 3, 23.5, 11, 60],  # cannedmilk
        [20.6, 17, .6, 6, 55.8, .2, 0, 0, 0],  # margarine
        [7.4, 448, 16.4, 19, 28.1, .8, 10.3, 4, 0],  # cheese
        [15.7, 661, 1, 48, 0, 9.6, 8.1, 471, 0],  # peanut butter
        [41.7, 0, 0, 0, .2, 0, .5, 5, 0],  # lard
        [2.2, 333, .2, 139, 169.2, 6.4, 50.8, 316, 525],  # liver
        [4.4, 249, .3, 37, 0, 18.2, 3.6, 79, 0],  # porkroast
        [5.8, 705, 6.8, 45, 3.5, 1, 4.9, 209, 0],  # salmon
        [2.4, 138, 3.7, 80, 69, 4.3, 5.8, 37, 862],  # greenbeans
        [2.6, 125, 4, 36, 7.2, 9, 4.5, 26, 5369],  # cabbage
        [5.8, 166, 3.8, 59, 16.6, 4.7, 5.9, 21, 1184],  # onions
        [14.3, 336, 1.8, 118, 6.7, 29.4, 7.1, 198, 2522],  # potatoes
        [1.1, 106, 0.0, 138, 918.4, 5.7, 13.8, 33, 2755],  # spinach
        [9.6, 138, 2.7, 54, 290.7, 8.4, 5.4, 83, 1912],  # sweet potatos
        [8.5, 87, 1.7, 173, 86.8, 1.2, 4.3, 55, 57],  # peaches
        [12.8, 99, 2.5, 154, 85.7, 3.9, 4.3, 65, 257],  # prunes
        [17.4, 1055, 3.7, 459, 5.1, 26.9, 38.2, 93, 0],  # limabeans
        [26.9, 1691, 11.4, 792, 0, 38.4, 24.6, 217, 0]] # navybeans

    daily_allowance = [
        3., 70., .8, 12., 5., 1.8, 2.7, 18., 75.]
    x,y = DietModel('Stinglers Diet Problem',
                    foods,
                    nutrients,
                    daily_allowance,
                    nutritive_value)

    # with DietModel('Stinglers Diet Problem',
    #                foods,
    #                nutrients,
    #                daily_allowance,
    #                nutritive_value) as M:
    #     M.solve()
    if True:
        print("Solution:")
        #x = M.dailyPurchase()
        for name, quantum in zip(foods, x):
            print("%20s : $%.4f " % (name, quantum))

        print("Nutrients:")
        #y = M.dailyNutrients()
        for name, quantum, unit, minval in zip(nutrients, y, nutrient_unit, daily_allowance):
            print("%20s : %7.2f %s (%.2f)" % (name, quantum, unit, minval))

djc1.py

Listing 16.9 djc1.py Click here to download.
##
#   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#   File:      djc1.py
#
#   Purpose: Demonstrates how to solve the problem with two disjunctions:
#
#      minimize    2x0 + x1 + 3x2 + x3
#      subject to   x0 + x1 + x2 + x3 >= -10
#                  (x0-2x1<=-1 and x2=x3=0) or (x2-3x3<=-2 and x1=x2=0)
#                  x0=2.5 or x1=2.5 or x2=2.5 or x3=2.5
##

from mosek.fusion import *
import mosek.fusion.pythonic
import sys

with Model('djc1') as M:
    # Create a variable of length 4
    x = M.variable('x', 4)

    # First disjunctive constraint
    M.disjunction( (x[0]-2*x[1]<=-1) & (x[2:4]==0) 
                  | 
                   (x[2]-3*x[3]<=-2) & (x[0:2]==0))
    
    # Second disjunctive constraint
    # Array of terms reading x_i = 2.5 for i = 0,1,2,3
    M.disjunction( [ x[i] == 2.5 for i in range(4) ])

    # The linear constraint
    M.constraint(Expr.sum(x) >= -10)

    # Objective
    M.objective(ObjectiveSense.Minimize, Expr.dot([2,1,3,1], x))

    # Useful for debugging
    M.writeTask("djc.ptf")          # Save to a readable file
    M.setLogHandler(sys.stdout)     # Enable log output

    # Solve 
    M.solve()

    if M.getPrimalSolutionStatus() == SolutionStatus.Optimal:
        print(f"Solution ={x.level()}")
    else:
        print("Another solution status")

duality.py

Listing 16.10 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.11 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 *
import mosek.fusion.pythonic

# Construct the model with parameters b, lambda1, lambda2
# and with prescribed matrix A
def initializeModel(m, n, A):
    M = Model()
    x = M.variable("x", n)

    # t >= |Ax-b|_2 where b is a parameter
    b = M.parameter("b", m)
    t = M.variable()
    M.constraint(Expr.vstack(t, A@x-b) == Domain.inQCone())

    # p_i >= |x_i|, i=1..n
    p = M.variable(n)
    M.constraint(Expr.hstack(p, x) == Domain.inQCone())

    # q >= |x|_2
    q = M.variable()
    M.constraint(Expr.vstack(q, x) == Domain.inQCone())

    # Objective, parametrized with lambda1, lambda2
    # t + lambda1*sum(p) + lambda2*q
    lambda1 = M.parameter("lambda1")
    lambda2 = M.parameter("lambda2")
    obj = t + lambda1 * Expr.sum(p) + lambda2 * q
    M.objective(ObjectiveSense.Minimize, obj)

    # Return the ready model
    return M

def smallExample():
    # Create a small example
    m, n = 4, 2
    A = [ [1.0, 2.0],
          [3.0, 4.0],
          [-2.0, -1.0],
          [-4.0, -3.0] ]
    M = initializeModel(m, n, A)

    # For convenience retrieve some elements of the model
    b = M.getParameter("b")
    lambda1 = M.getParameter("lambda1")
    lambda2 = M.getParameter("lambda2")
    x = M.getVariable("x")

    # First solve
    b.setValue([0.1, 1.2, -1.1, 3.0])
    lambda1.setValue(0.1)
    lambda2.setValue(0.01)

    M.solve()
    print("Objective {0}, solution x = {1}".format(M.primalObjValue(), x.level()))
    
    # Increase lambda1
    lambda1.setValue(0.5)
    
    M.solve()
    print("Objective {0}, solution x = {1}".format(M.primalObjValue(), x.level()))

    # Now change the data completely
    b.setValue([1.0, 1.0, 1.0, 1.0])
    lambda1.setValue(0.0)
    lambda2.setValue(0.0)
    
    M.solve()
    print("Objective {0}, solution x = {1}".format(M.primalObjValue(), x.level()))

    # And increase lamda2
    lambda2.setValue(1.4145)
    
    M.solve()
    print("Objective {0}, solution x = {1}".format(M.primalObjValue(), x.level()))

    M.dispose()

smallExample()

facility_location.py

Listing 16.12 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.13 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 mosek.fusion.pythonic
import sys

# Models log(sum(exp(Ax+b))) <= 0.
# Each row of [A b] describes one of the exp-terms
def logsumexp(M, A, x, b):    
    k = int(A.shape[0])
    u = M.variable(k)
    M.constraint(Expr.sum(u) == 1.0)
    M.constraint(Expr.hstack(u,
                             Expr.constTerm(k, 1.0),
                             A @ x + b), Domain.inPExpCone())

# maximize     h*w*d
# subjecto to  2*(h*w + h*d) <= Awall
#              w*d <= Afloor
#              alpha <= h/w <= beta
#              gamma <= d/w <= delta
#
# Variable substitutions:  h = exp(x), w = exp(y), d = exp(z).
#
# maximize     x+y+z
# subject      log( exp(x+y+log(2/Awall)) + exp(x+z+log(2/Awall)) ) <= 0
#                              y+z <= log(Afloor)
#              log( alpha ) <= x-y <= log( beta )
#              log( gamma ) <= z-y <= log( delta )
def max_volume_box(Aw, Af, alpha, beta, gamma, delta):
    with Model('max_vol_box') as M:
        xyz = M.variable(3)
        M.objective('Objective', ObjectiveSense.Maximize, Expr.sum(xyz))
        
        logsumexp(M, array([[1,1,0],[1,0,1]]), xyz, array([log(2.0/Aw), log(2.0/Aw)]))
        
        M.constraint(Expr.dot([0, 1, 1], xyz) <= log(Af))
        M.constraint(Expr.dot([1,-1, 0], xyz) == Domain.inRange(log(alpha),log(beta)))
        M.constraint(Expr.dot([0,-1, 1], xyz) == Domain.inRange(log(gamma),log(delta)))
        
        M.setLogHandler(sys.stdout)
        M.solve()
        
        return exp(xyz.level())

Aw, Af, alpha, beta, gamma, delta = 200.0, 50.0, 2.0, 10.0, 2.0, 10.0
h,w,d = max_volume_box(Aw, Af, alpha, beta, gamma, delta) 
print("h={0:.3f}, w={1:.3f}, d={2:.3f}".format(h, w, d))

lo1.py

Listing 16.14 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 *
import mosek.fusion.pythonic

def main(args):
    A = [[3.0, 1.0, 2.0, 0.0],
         [2.0, 1.0, 3.0, 1.0],
         [0.0, 2.0, 0.0, 3.0]]
    c = [3.0, 1.0, 5.0, 1.0]

    # Create a model with the name 'lo1'
    with Model("lo1") as M:

        # Create variable 'x' of length 4
        x = M.variable("x", 4, Domain.greaterThan(0.0))

        # Create constraints
        M.constraint(x[1] <= 10)
        M.constraint("c1", x.T @ A[0] == 30.0)
        M.constraint("c2", x.T @ A[1] >= 15.0)
        M.constraint("c3", x.T @ A[2] <= 25.0)

        # Set the objective function to (c^t * x)
        M.objective("obj", ObjectiveSense.Maximize, x.T @ c)

        # Solve the problem
        M.solve()

        # Get the solution values
        sol = x.level()
        print('\n'.join(["x[%d] = %f" % (i, sol[i]) for i in range(4)]))

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

logistic.py

Listing 16.15 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 mosek.fusion.pythonic
import numpy as np 
import sys, itertools

# t >= log( 1 + exp(u) ) coordinatewise
def softplus(M, t, u):
    n = t.getShape()[0]
    z1 = M.variable(n)
    z2 = M.variable(n)
    M.set(z1 + z2 == 1,
          Expr.hstack(z1, Expr.constTerm(n, 1.0), u-t) == Domain.inPExpCone(),
          Expr.hstack(z2, Expr.constTerm(n, 1.0), -t)  == Domain.inPExpCone())

# Model logistic regression (regularized with full 2-norm of theta)
# X - n x d matrix of data points
# y - length n vector classifying training points
# lamb - regularization parameter
def logisticRegression(X, y, lamb=1.0):
    n, d = int(X.shape[0]), int(X.shape[1])         # num samples, dimension
    M = Model()
    theta = M.variable(d)
    t     = M.variable(n)
    reg   = M.variable()

    M.objective(ObjectiveSense.Minimize, Expr.sum(t) + lamb * reg)
    M.constraint(Var.vstack(reg, theta), Domain.inQCone())

    signs = list(map(lambda y: -1.0 if y==1 else 1.0, y))
    softplus(M, t, Expr.mulElm(X @ theta, signs))

    return M, theta

# Map the 2d data through all monomials of degree <= d
def mapFeature(p,d):
    return np.array([(p[0]**a)*(p[1]**b) for a,b in itertools.product(range(d+1), range(d+1)) if a+b<=d])
def mapFeatures(x,d):
    return np.array([mapFeature(p,d) for p in x])

# Load the file and map using degree d monomials
# The file format is
# x_1, x_2, y
# for all training examples
def loaddata(filename):
    # Read coordinates and y values
    x, y = [], []
    with open(filename, "r") as f:
        for l in f.readlines():
            num = l.split(',')
            x.append([float(num[0]), float(num[1])])
            y.append(int(num[2]))
    return np.array(x), np.array(y)

# Plot some 2d results
def plot2d(x, y, d, theta):
    import matplotlib
    import matplotlib.pyplot as plt

    pos = np.where(y==1)
    neg = np.where(y==0)
    plt.scatter(x[pos,0], x[pos,1], marker='o', color='b')
    plt.scatter(x[neg,0], x[neg,1], marker='x', color='r')

    u = np.linspace(-1, 1, 50)
    v = np.linspace(-1, 1, 50)
    z = np.zeros(shape=(len(u), len(v)))
    for i in range(len(u)):
        for j in range(len(v)):
            z[i,j] = np.dot(mapFeature([u[i],v[j]], d), theta)
    plt.contour(u, v, z.T, [0])

    plt.show()

###############################################################################

# Example from documentation is contained in
# https://datascienceplus.com/wp-content/uploads/2017/02/ex2data2.txt
'''
for lamb in [1e-0,1e-2,1e-4]:
    x, y = loaddata("ex2data2.txt")
    X = mapFeatures(x, d=6)
    M, theta = logisticRegression(X, y, lamb)
    M.solve()    
    plot2d(x, y, 6, theta.level())
'''

# Example 2: discover a circle
x, y = [], []
for x1,x2 in itertools.product(np.linspace(-1, 1, 30),np.linspace(-1, 1, 30)):
    x.append([x1,x2])
    y.append(0 if x1**2+x2**2<0.69 else 1)
x, y = np.array(x), np.array(y)
X = mapFeatures(x, d=2)
M, theta = logisticRegression(X, y, 0.1)
M.setLogHandler(sys.stdout)

M.solve()
try:
    print(theta.level())
except:
    print("Could not get solution")
    sys.exit(1)

lownerjohn_ellipsoid.py

Listing 16.16 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 *
import mosek.fusion.pythonic

'''
 Purpose: Models the hypograph of the n-th power of the
 determinant of a positive definite matrix. See [1,2] for more details.

   The convex set (a hypograph)

   C = { (X, t) \\in S^n_+ x R |  t <= det(X)^{1/n} },

   can be modeled as the intersection of a semidefinite cone

   [ X, Z; Z^T Diag(Z) ] >= 0

   and a geometric mean bound

   t <= (Z11*Z22*...*Znn)^{1/n} 
'''
def det_rootn(M, t, n):
    # Setup variables
    Y = M.variable(Domain.inPSDCone(2 * n))

    # Setup Y = [X, Z; Z^T , diag(Z)]
    X   = Y[0:n, 0:n]
    Z   = Y[0:n, n:2*n]
    DZ  = Y[n:2*n, n:2*n]

    # Z is lower-triangular
    M.constraint(Z.pick([[i,j] for i in range(n) for j in range(i+1,n)]), Domain.equalsTo(0.0))
    # DZ = Diag(Z)
    M.constraint(DZ == Expr.mulElm(Z, Matrix.eye(n)))

    # (Z11*Z22*...*Znn) >= t^n
    M.constraint(Expr.vstack(DZ.diag(), t), Domain.inPGeoMeanCone())

    # Return an n x n PSD variable which satisfies t <= det(X)^(1/n)
    return X

'''
  The inner ellipsoidal approximation to a polytope

     S = { x \\in R^n | Ax < b }.

  maximizes the volume of the inscribed ellipsoid,

     { x | x = C*u + d, || u ||_2 <= 1 }.

  The volume is proportional to det(C)^(1/n), so the
  problem can be solved as

    maximize         t
    subject to       t       <= det(C)^(1/n)
                || C*ai ||_2 <= bi - ai^T * d,  i=1,...,m
                C is PSD

  which is equivalent to a mixed conic quadratic and semidefinite
  programming problem.
'''
def lownerjohn_inner(A, b):
    with Model("lownerjohn_inner") as M:
        M.setLogHandler(sys.stdout)
        m, n = len(A), len(A[0])

        # Setup variables
        t = M.variable("t", 1, Domain.greaterThan(0.0))
        C = det_rootn(M, t, n)
        d = M.variable("d", n, Domain.unbounded())

        # (b-Ad, AC) generate cones
        M.constraint("qc", Expr.hstack(b - A @ d, A @ C), Domain.inQCone())

        # Objective: Maximize t
        M.objective(ObjectiveSense.Maximize, t)

        M.writeTask("lj-inner.ptf")
        M.solve()

        C, d = C.level(), d.level()
        return ([C[i:i + n] for i in range(0, n * n, n)], d)

'''
  The outer ellipsoidal approximation to a polytope given
  as the convex hull of a set of points

    S = conv{ x1, x2, ... , xm }

  minimizes the volume of the enclosing ellipsoid,

    { x | || P*x-c ||_2 <= 1 }

  The volume is proportional to det(P)^{-1/n}, so the problem can
  be solved as

    maximize         t
    subject to       t       <= det(P)^(1/n)
                || P*xi - c ||_2 <= 1,  i=1,...,m
                P is PSD.
'''
def lownerjohn_outer(x):
    with Model("lownerjohn_outer") as M:
        M.setLogHandler(sys.stdout)
        m, n = len(x), len(x[0])

        # Setup variables
        t = M.variable("t", 1, Domain.greaterThan(0.0))
        P = det_rootn(M, t, n)
        c = M.variable("c", Domain.unbounded().withShape(1,n))

        # (1, Px-c) in cone
        M.constraint("qc", Expr.hstack(Expr.ones(m), x @ P - Expr.repeat(c,m,0)), Domain.inQCone())

        # Objective: Maximize t
        M.objective(ObjectiveSense.Maximize, t)
        M.writeTask("lj-outer.ptf")
        M.solve()

        P, c = P.level(), c.level()
        return ([P[i:i + n] for i in range(0, n * n, n)], c)

##########################################################################

if __name__ == '__main__':
    #Vertices of a polygon in 2D
    p = [[0., 0.], [1., 3.], [5.5, 4.5], [7., 4.], [7., 1.], [3., -2.]]
    nVerts = len(p)

    #The hyperplane representation of the same polytope
    A = [[-p[i][1] + p[i - 1][1], p[i][0] - p[i - 1][0]]
         for i in range(len(p))]
    b = [A[i][0] * p[i][0] + A[i][1] * p[i][1] for i in range(len(p))]

    Po, co = lownerjohn_outer(p)
    Ci, di = lownerjohn_inner(A, b)


    #Visualization
    try:
        import numpy as np
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        import matplotlib.patches as patches

        #Polygon
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.add_patch(patches.Polygon(p, fill=False, color="red"))
        #The inner ellipse
        theta = np.linspace(0, 2 * np.pi, 100)
        x = Ci[0][0] * np.cos(theta) + Ci[0][1] * np.sin(theta) + di[0]
        y = Ci[1][0] * np.cos(theta) + Ci[1][1] * np.sin(theta) + di[1]
        ax.plot(x, y)
        #The outer ellipse
        x, y = np.meshgrid(np.arange(-1.0, 8.0, 0.025),
                           np.arange(-3.0, 6.5, 0.025))
        ax.contour(x, y,
                   (Po[0][0] * x + Po[0][1] * y - co[0])**2 + (Po[1][0] * x + Po[1][1] * y - co[1])**2, [1])
        ax.autoscale_view()
        ax.xaxis.set_visible(False)
        ax.yaxis.set_visible(False)
        fig.savefig('ellipsoid.png')
    except:
        print("Inner:")
        print("  C = ", Ci)
        print("  d = ", di)
        print("Outer:")
        print("  P = ", Po)
        print("  c = ", co)

lpt.py

Listing 16.17 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.18 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 *
import mosek.fusion.pythonic

with Model('mico1') as M:

    x = M.variable(Domain.integral(Domain.unbounded()))
    y = M.variable(Domain.integral(Domain.unbounded()))
    t = M.variable()

    M.constraint(Expr.vstack(t, x, y), Domain.inQCone())
    M.constraint(Expr.vstack(x - 3.8, 1, y), Domain.inPExpCone())

    M.objective(ObjectiveSense.Minimize, t)

    M.setLogHandler(sys.stdout)
    M.solve()

    print('Solution: x = {0}, y = {1}'.format(x.level()[0], y.level()[0]))

milo1.py

Listing 16.19 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.20 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)

        # Request constructing the solution from integer variable values
        M.setSolverParam("mioConstructSol", "on")

        # Solve the problem
        M.setLogHandler(sys.stdout)
        M.solve()

        # Get the solution values
        ss = M.getPrimalSolutionStatus()
        print("Solution status: {0}".format(ss))
        sol = x.level()
        print('x = {0}'.format(sol))

        # Was the initial solution used?
        constr = M.getSolverIntInfo("mioConstructSolution")
        constrVal = M.getSolverDoubleInfo("mioConstructSolutionObj")
        print("Construct solution utilization: {0}\nConstruct solution objective: {1:.3f}\n".format(constr, constrVal))                    

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

nearestcorr.py

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 *
import mosek.fusion.pythonic
from mosek import LinAlg

"""
 Assuming that e is an NxN expression, return the lower triangular part as a vector.
"""
def vec(e):
    N = e.getShape()[0]

    msubi = range(N * (N + 1) // 2)
    msubj = [i * N + j for i in range(N) for j in range(i + 1)]
    mcof  = [2.0**0.5 if i !=
             j else 1.0 for i in range(N) for j in range(i + 1)]

    S = Matrix.sparse(N * (N + 1) // 2, N * N, msubi, msubj, mcof)
    return S @ Expr.flatten(e)

def nearestcorr(A):
    N = A.numRows()

    # Create a model
    with Model("NearestCorrelation") as M:
        # Setting up the variables
        X = M.variable("X", Domain.inPSDCone(N))
        t = M.variable("t", 1, Domain.unbounded())

        # (t, vec (A-X)) \in Q
        M.constraint("C1", Expr.vstack(t, vec(A-X)), Domain.inQCone())

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

        # Objective: Minimize t
        M.objective(ObjectiveSense.Minimize, t)
        M.writeTask('nearcor.ptf')
        M.solve()

        return X.level(), t.level()

def nearestcorr_nucnorm(A, gammas):
    N = A.numRows()
    with Model("NucNorm") as M:
        # Setup variables
        t = M.variable("t", 1, Domain.unbounded())
        X = M.variable("X", Domain.inPSDCone(N))
        w = M.variable("w", N, Domain.greaterThan(0.0))

        # D = diag(w)
        D = Expr.mulElm(Matrix.eye(N), Var.repeat(w, 1, N))
        # (t, vec (X + D - A)) in Q
        M.constraint(Expr.vstack(t, vec(X+D-A)), Domain.inQCone())

        result = []
        for g in gammas:
            # Objective: Minimize t + gamma*Tr(X)
            M.objective(ObjectiveSense.Minimize, t + g * Expr.sum(X.diag()))
            M.solve()

            # Find eigenvalues of X and compute its rank
            d = [0.0] * int(N)
            LinAlg.syeig(mosek.uplo.lo, N, X.level(), d)
            result.append(
                (g, t.level()[0], sum([d[i] > 1e-6 for i in range(N)]), X.level()))

        return result

if __name__ == '__main__':
    N = 5
    A = Matrix.dense(N, N, [0.0, 0.5, -0.1, -0.2, 0.5,
                            0.5, 1.25, -0.05, -0.1, 0.25,
                            -0.1, -0.05, 0.51, 0.02, -0.05,
                            -0.2, -0.1, 0.02, 0.54, -0.1,
                            0.5, 0.25, -0.05, -0.1, 1.25])

    gammas = [0.1 * i for i in range(11)]

    X, t = nearestcorr(A)

    print("--- Nearest Correlation ---")
    print("X = ")
    print(X.reshape((N, N)))
    print("t = ", t)

    print("--- Nearest Correlation with Nuclear Norm---")
    for g, res, rank, X in nearestcorr_nucnorm(A, gammas):
        print("gamma=%f, res=%e, rank=%d" % (g, res, rank))

opt_server_sync.py

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) < 2:
    print("Missing argument, syntax is:")
    print("  python opt_server_sync.py addr [certpath]")
    sys.exit(1)

serveraddr = sys.argv[1]
tlscert   = None if len(sys.argv)<3 else sys.argv[2]

with Model('testOptServer') as M:
    # Setup a simple test problem
    x = M.variable('x', 3, Domain.greaterThan(0.0))
    M.constraint("lc", Expr.dot([1.0, 1.0, 2.0], x), Domain.equalsTo(1.0))
    M.objective("obj", ObjectiveSense.Minimize, Expr.sum(x))

    # Attach log handler
    M.setLogHandler(sys.stdout)

    # Set OptServer URL
    M.optserverHost(serveraddr)

    # Path to certificate, if any
    M.setSolverParam("remoteTlsCertPath", tlscert)

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

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

parallel.py

Listing 16.23 parallel.py Click here to download.
#
#  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File:      parallel.py
#
#  Purpose: Demonstrates parallel optimization using solveBatch()
#
from mosek.fusion import *
 
def makeToyParameterizedModel():
  M = Model()
  x = M.variable("x",3)
  p = M.parameter("p")
  M.objective(ObjectiveSense.Maximize, Expr.sum(x))
  M.constraint(Expr.vstack(p,x), Domain.inQCone())
  return M

# Example of how to use Model.solveBatch()
def main():
  # Choose some sample parameters
  n = 10                 # Number of models to optimize
  threadpoolsize = 4     # Total number of threads available
  threadspermodel = 1    # Number of threads per each model

  # Create a toy model for this example
  M = makeToyParameterizedModel()

  # Set up n copies of the model with different data
  models = [M.clone() for _ in range(n)]
  for i in range(n):
    models[i].getParameter("p").setValue(i+1)
    # We can set the number of threads individually per model
    models[i].setSolverParam("numThreads", threadspermodel)

  # Solve all models in parallel
  status = Model.solveBatch(False,         # No race
                            -1.0,          # No time limit
                            threadpoolsize,
                            models)        # Array of Models to solve

  # Access information about each model
  for i in range(n):
    if status[i] == SolverStatus.OK:
      print("Model {}: Status {}, Solution Status {}, Objective {:.3f}, Time {:.3f}".format(
          i, 
          status[i],
          models[i].getPrimalSolutionStatus(),
          models[i].primalObjValue(),
          models[i].getSolverDoubleInfo("optimizerTime")))
    else:
      print("Model {}: not solved".format(i))

if __name__ == "__main__":
    main()

parameters.py

Listing 16.24 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))

pinfeas.py

Listing 16.25 pinfeas.py Click here to download.
##
#   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#   File:      pinfeas.py
#
#   Purpose: Demonstrates how to fetch a primal infeasibility certificate
#            for a linear problem
#
##

from mosek.fusion import *
import sys

# Analyzes and prints infeasibility certificate for a single object,
# which can be a variable or constraint
def analyzeCertificate(name,          # name of the analyzed object
                       size,          # size of the object
                       duals,         # double[], actual dual values
                       eps):          # tolerance determining when a dual value is considered important
    for i in range(0, size):
        if abs(duals[i]) > eps:
            print(f"{name}[{i}], dual = {duals[i]}")

# Construct the sample model from the example in the manual
sMat = Matrix.sparse(3, 7, [0,0,1,1,2,2,2],
                           [0,1,2,3,4,5,6],
                           [1.0]*7)
sBound = [200, 1000, 1000]
dMat = Matrix.sparse(4, 7, [0,0,1,2,2,3,3],
                           [0,4,1,2,5,3,6],
                           [1.0]*7)
dBound = [1100, 200, 500, 500]
c = [1, 2, 5, 2, 1, 2, 1]

M = Model("pinfeas")
x = M.variable("x", 7, Domain.greaterThan(0))
s = M.constraint("s", Expr.mul(sMat, x), Domain.lessThan(sBound))
d = M.constraint("d", Expr.mul(dMat, x), Domain.equalsTo(dBound))
M.objective(ObjectiveSense.Minimize, Expr.dot(c,x))

# Set useful debugging tools and solve the model
M.setLogHandler(sys.stdout)
M.writeTask("pinfeas.ptf")
M.solve()

# Check problem status
if M.getProblemStatus() == ProblemStatus.PrimalInfeasible:
    # Set the tolerance at which we consider a dual value as essential
    eps = 1e-7

    # We want to retrieve infeasibility certificates
    M.acceptedSolutionStatus(AccSolutionStatus.Certificate)

    # Go through variable bounds
    print("Variable bounds important for infeasibility: ")
    analyzeCertificate(name = "x", size = x.getSize(), duals = x.dual(), eps = eps)

    # Go through constraint bounds
    print("Constraint bounds important for infeasibility: ")
    analyzeCertificate(name = "s", size = s.getSize(), duals = s.dual(), eps = eps)
    analyzeCertificate(name = "d", size = d.getSize(), duals = d.dual(), eps = eps)

else:
    print("The problem is not primal infeasible, no certificate to show")

portfolio_1_basic.py

Listing 16.26 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 mosek.fusion.pythonic
import numpy as np

"""
Purpose:
    Computes the optimal portfolio for a given risk 
 
Input:
    n: Number of assets
    mu: An n dimensional vector of expected returns
    GT: A matrix with n columns so (GT')*GT  = covariance matrix
    x0: Initial holdings 
    w: Initial cash holding
    gamma: Maximum risk (=std. dev) accepted
 
Output:
    Optimal expected return and the optimal portfolio     
""" 
def BasicMarkowitz(n,mu,GT,x0,w,gamma):
    
    with  Model("Basic Markowitz") as M:

        # Redirect log output from the solver to stdout for debugging. 
        # if uncommented.
        # M.setLogHandler(sys.stdout) 
        
        # Defines the variables (holdings). Shortselling is not allowed.
        x = M.variable("x", n, Domain.greaterThan(0.0))
        
        #  Maximize expected return
        M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)
        
        # The amount invested  must be identical to initial wealth
        M.constraint('budget', Expr.sum(x) == w+sum(x0))
        
        # Imposes a bound on the risk
        M.constraint('risk', Expr.vstack(gamma, GT @ x) == Domain.inQCone())
        #M.constraint('risk', Expr.vstack(gamma, 0.5, Expr.mul(GT, x)), Domain.inRotatedQCone())

        # Solves the model.
        M.solve()

        # Check if the solution is an optimal point
        solsta = M.getPrimalSolutionStatus()
        if (solsta != SolutionStatus.Optimal):
            # See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
            raise Exception(f"Unexpected solution status: {solsta}")

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


if __name__ == '__main__':    

    n      = 8
    w      = 59
    mu     = [0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628]
    x0     = [8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0]
    gammas = [36]
    GT     = [
        [0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638],
        [0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506],
        [0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914],
        [0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742],
        [0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ]
    ]


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

portfolio_2_frontier.py

Listing 16.27 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 mosek.fusion.pythonic
import numpy as np
import sys

"""
Purpose:
    Computes several portfolios on the optimal portfolios by

        for alpha in alphas:
            maximize   expected return - alpha * variance
            subject to the constraints

Input:
    n: Number of assets
    mu: An n dimensional vector of expected returns
    GT: A matrix with n columns so (GT')*GT  = covariance matrix
    x0: Initial holdings
    w: Initial cash holding
    alphas: List of the alphas

Output:
    The efficient frontier as list of tuples (alpha, expected return, variance)
"""
def EfficientFrontier(n,mu,GT,x0,w,alphas):

    with Model("Efficient frontier") as M:
        frontier = []

        # Defines the variables (holdings). Shortselling is not allowed.
        x = M.variable("x", n, Domain.greaterThan(0.0)) # Portfolio variables
        s = M.variable("s", 1, Domain.unbounded())      # Variance variable

        # Total budget constraint
        M.constraint('budget', Expr.sum(x) == w+sum(x0))

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

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

            M.solve()

            # Check if the solution is an optimal point
            solsta = M.getPrimalSolutionStatus()
            if (solsta != SolutionStatus.Optimal):
                # See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
                raise Exception(f"Unexpected solution status: {solsta}")

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

        return frontier


if __name__ == '__main__':

    n = 8
    w = 1.0   
    mu = [0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379]
    x0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    GT = [
        [0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638],
        [0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506],
        [0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914],
        [0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742],
        [0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ]
    ]

    # Some predefined alphas are chosen
    alphas = [0.0, 0.01, 0.1, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 10.0]
    frontier= EfficientFrontier(n,mu,GT,x0,w,alphas)
    print("\n-----------------------------------------------------------------------------------")
    print('Efficient frontier')
    print("-----------------------------------------------------------------------------------\n")
    print('%-12s  %-12s  %-12s' % ('alpha','return','risk (std. dev.)'))
    for i in frontier:
        print('%-12.4f  %-12.4e  %-12.4e' % (i[0],i[1],np.sqrt(i[2])))

portfolio_3_impact.py

Listing 16.28 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 mosek.fusion.pythonic
import numpy as np

"""
    Description:
        Extends the basic Markowitz model with a market cost term.

    Input:
        n: Number of assets
        mu: An n dimensional vector of expected returns
        GT: A matrix with n columns so (GT')*GT  = covariance matrix
        x0: Initial holdings 
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted
        m: It is assumed that  market impact cost for the j'th asset is
           m_j|x_j-x0_j|^3/2

    Output:
       Optimal expected return and the optimal portfolio     

"""
def MarkowitzWithMarketImpact(n,mu,GT,x0,w,gamma,m):
    with  Model("Markowitz portfolio with market impact") as M:

        #M.setLogHandler(sys.stdout) 
    
        # Defines the variables. No shortselling is allowed.
        x = M.variable("x", n, Domain.greaterThan(0.0))
        
        # Variables computing market impact 
        t = M.variable("t", n, Domain.unbounded())

        # Maximize expected return
        M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)

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

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

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

        M.solve()

        # Check if the solution is an optimal point
        solsta = M.getPrimalSolutionStatus()
        if (solsta != SolutionStatus.Optimal):
            # See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
            raise Exception(f"Unexpected solution status: {solsta}")

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


if __name__ == '__main__':    

    n = 8
    w = 1.0   
    mu = [0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379]
    x0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    GT = [
        [0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638],
        [0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506],
        [0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914],
        [0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742],
        [0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ]
    ]
                  
    # Somewhat arbitrary choice of m
    gamma = 0.36
    m = n * [0.01]
    xsol, tsol = MarkowitzWithMarketImpact(n,mu,GT,x0,w,gamma,m)
    print("\n-----------------------------------------------------------------------------------")
    print('Markowitz portfolio optimization with market impact cost')
    print("-----------------------------------------------------------------------------------\n")
    print('Expected return: %.4e Std. deviation: %.4e Market impact cost: %.4e' % \
          (np.dot(mu,xsol),gamma,np.dot(m,tsol)))
    print('Optimal portfolio: [%.4e, %.4e, %.4e, %.4e, %.4e, %.4e, %.4e, %.4e]' % tuple(xsol))

portfolio_4_transcost.py

Listing 16.29 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 mosek.fusion.pythonic
import numpy as np

"""
    Description:
        Extends the basic Markowitz model with a market cost term.

    Input:
        n: Number of assets
        mu: An n dimensional vector of expected returns
        GT: A matrix with n columns so (GT')*GT  = covariance matrix
        x0: Initial holdings 
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted
        f: If asset j is traded then a fixed cost f_j must be paid
        g: If asset j is traded then a cost g_j must be paid for each unit traded

    Output:
       Optimal expected return and the optimal portfolio     

"""
def MarkowitzWithTransactionsCost(n,mu,GT,x0,w,gamma,f,g):
    # Upper bound on the traded amount
    w0 = w+sum(x0)
    u = n*[w0]

    with Model("Markowitz portfolio with transaction costs") as M:
        #M.setLogHandler(sys.stdout)

        # Defines the variables. No shortselling is allowed.
        x = M.variable("x", n, Domain.greaterThan(0.0))

        # Additional "helper" variables 
        z = M.variable("z", n, Domain.unbounded())   
        # Binary variables
        y = M.variable("y", n, Domain.binary())

        #  Maximize expected return
        M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)

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

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

        # z >= |x-x0| 
        M.constraint('buy',  z >= x-x0)
        M.constraint('sell', z >= x0-x)
        # Alternatively, formulate the two constraints as
        #M.constraint('trade', Expr.hstack(z, x-x0), Domain.inQcone())

        # Constraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
        M.constraint('y_on_off', z <= Expr.mulElm(u, y))

        # Integer optimization problems can be very hard to solve so limiting the 
        # maximum amount of time is a valuable safe guard
        M.setSolverParam('mioMaxTime', 180.0) 
        M.solve()

        # Check if the solution is an optimal point
        solsta = M.getPrimalSolutionStatus()
        if (solsta != SolutionStatus.Optimal):
            # See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
            raise Exception(f"Unexpected solution status: {solsta}")

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

if __name__ == '__main__':    

    n = 8
    w = 1.0   
    mu = [0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379]
    x0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    GT = [
        [0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638],
        [0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506],
        [0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914],
        [0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742],
        [0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ]
    ]

    f = n*[0.01]
    g = n*[0.001]
    gamma = 0.36
    xsol, ysol, zsol = MarkowitzWithTransactionsCost(n,mu,GT,x0,w,gamma,f,g)
    print("\n-----------------------------------------------------------------------------------")
    print('Markowitz portfolio optimization with transactions cost')
    print("-----------------------------------------------------------------------------------\n")
    print('Expected return: %.4e Std. deviation: %.4e Transactions cost: %.4e' % \
          (np.dot(mu,xsol),gamma, np.dot(f,ysol)+np.dot(g,zsol)))

portfolio_5_card.py

Listing 16.30 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 mosek.fusion.pythonic
import numpy as np 

"""
    Description:
        Extends the basic Markowitz model with cardinality constraints.

    Input:
        n: Number of assets
        mu: An n dimensional vector of expected returns
        GT: A matrix with n columns so (GT')*GT  = covariance matrix
        x0: Initial holdings 
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted
        k: Maximum number of assets on which we allow to change position.

    Output:
       Optimal expected return and the optimal portfolio     

"""
def MarkowitzWithCardinality(n,mu,GT,x0,w,gamma,KValues):
    # Upper bound on the traded amount
    w0 = w+sum(x0)
    u = n*[w0]

    with Model("Markowitz portfolio with cardinality bound") as M:
        #M.setLogHandler(sys.stdout)

        # Defines the variables. No shortselling is allowed.
        x = M.variable("x", n, Domain.greaterThan(0.0))

        # Additional "helper" variables 
        z = M.variable("z", n, Domain.unbounded())   
        # Binary variables  - do we change position in assets
        y = M.variable("y", n, Domain.binary())

        #  Maximize expected return
        M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)

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

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

        # z >= |x-x0| 
        M.constraint('buy',  z >= x-x0)
        M.constraint('sell', z >= x0-x)

        # Constraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
        M.constraint('y_on_off', z <= Expr.mulElm(u,y))

        # At most K assets change position
        cardMax = M.parameter()
        M.constraint('cardinality', Expr.sum(y)<= cardMax)

        # Integer optimization problems can be very hard to solve so limiting the 
        # maximum amount of time is a valuable safe guard
        M.setSolverParam('mioMaxTime', 180.0) 

        # Solve multiple instances by varying the parameter K
        results = []
        for K in KValues:
            cardMax.setValue(K)
            M.solve()
            
            # Check if the solution is an optimal point
            solsta = M.getPrimalSolutionStatus()
            if (solsta != SolutionStatus.Optimal):
                # See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
                raise Exception(f"Unexpected solution status: {solsta}")

            results.append(x.level())

        return results


if __name__ == '__main__':    

    n = 8
    w = 1.0   
    mu = [0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379]
    x0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    GT = [
        [0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638],
        [0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506],
        [0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914],
        [0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742],
        [0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327],
        [0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ]
    ]
    gamma  = 0.25

    xsol = MarkowitzWithCardinality(n,mu,GT,x0,w,gamma,range(1,n+1)) 
    print("\n-----------------------------------------------------------------------------------")
    print('Markowitz portfolio optimization with cardinality constraints')
    print("-----------------------------------------------------------------------------------\n")
    for K in range(0,n):
        print('Bound: {0}   Expected return: {1:.4f}  Solution {2}'.format(K+1, np.dot(mu, xsol[K]), xsol[K]))

portfolio_6_factor.py

Listing 16.31 portfolio_6_factor.py Click here to download.
##
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File:      portfolio_6_factor.py
#
#  Description :  Implements a basic portfolio optimization model
#                 with factor structured covariance matrix.
##

from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np 

"""
    Description:
        Extends the basic Markowitz model with factor structure in the covariance matrix.

    Input:
        n: Number of securities
        mu: An n dimensional vector of expected returns
        G_factor_T: The factor (dense) part of the factorized risk
        theta: specific risk vector
        x0: Initial holdings 
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted

    Output:
       Optimal expected return and the optimal portfolio     

"""
def FactorModelMarkowitz(n, mu, G_factor_T, S_theta, x0, w, gamma):
    
    with Model("Factor model Markowitz") as M:

        # Variables 
        # The variable x is the fraction of holdings in each security. 
        # It is restricted to be positive, which imposes the constraint of no short-selling.   
        x = M.variable("x", n, Domain.greaterThan(0.0))

        # Objective (quadratic utility version)
        M.objective('obj', ObjectiveSense.Maximize, x.T @ mu)
       
        # Budget constraint
        M.constraint('budget', Expr.sum(x) == w + sum(x0))

        # Conic constraint for the portfolio std. dev
        M.constraint('risk', Expr.vstack([Expr.constTerm(gamma),
                                          G_factor_T @ x, 
                                          Expr.mulElm(np.sqrt(theta), x)]), Domain.inQCone())

        # Solve optimization
        M.solve()

        # Check if the solution is an optimal point
        solsta = M.getPrimalSolutionStatus()
        if (solsta != SolutionStatus.Optimal):
            # See https://docs.mosek.com/latest/pythonfusion/accessing-solution.html about handling solution statuses.
            raise Exception(f"Unexpected solution status: {solsta}")
        
        return mu @ x.level(), x.level()

if __name__ == '__main__':

    n = 8
    w = 1.0   
    mu = [0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379]
    x0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    B = np.array([
        [0.4256, 0.1869],
        [0.2413, 0.3877],
        [0.2235, 0.3697],
        [0.1503, 0.4612],
        [1.5325, -0.2633],
        [1.2741, -0.2613],
        [0.6939, 0.2372],
        [0.5425, 0.2116]
    ])
    S_F = np.array([
        [0.0620, 0.0577],
        [0.0577, 0.0908]
    ])
    theta = np.array([0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459])
    P           = np.linalg.cholesky(S_F)
    G_factor    = B @ P
    G_factor_T  = G_factor.T

    gammas = [0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48] 

    print("\n-----------------------------------------------------------------------------------")
    print('Markowitz portfolio optimization with factor model')
    print("-----------------------------------------------------------------------------------\n")
    for gamma in gammas:
        er, x = FactorModelMarkowitz(n, mu, G_factor_T, theta, x0, w, gamma)
        print('Expected return: %.4e Std. deviation: %.4e' % (er, gamma))
        np.set_printoptions(precision=4)
        print(f'Optimal portfolio: {x}')

pow1.py

Listing 16.32 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.33 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 mosek.fusion.pythonic

import random 

def primal_svm(m,n,X,y,CC):

    print("Number of data    : %d"%m)
    print("Number of features: %d"%n)

    with Model() as M:

        w =  M.variable('w' , n, Domain.unbounded())
        t =  M.variable('t' , 1, Domain.unbounded())
        b =  M.variable('b' , 1, Domain.unbounded())
        xi = M.variable('xi', m, Domain.greaterThan(0.))

        M.constraint( Expr.mulElm( y, X @ w - Var.repeat(b,m) ) >= 1 - xi )

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

        print ('   c   |    b      |     w ')

        for C in CC:
            M.objective( ObjectiveSense.Minimize, t + C * Expr.sum(xi) )
            M.solve()

            try:
                cb = '{0:6} | {1:8f} | '.format(C,b.level()[0]) 
                wstar =' '.join([ '{0:8f}'.format(wi) for wi in w.level()]) 
                print (cb+wstar)
            except:
                pass;             

if __name__ == '__main__':

    CC=[ 500.0*i for i in range(10)]

    m  = 50
    n  = 3
    seed= 0

    random.seed(seed)
    nump= random.randint(0,50)
    numm= m - nump

    y = [  1. for i in range(nump)] + \
        [ -1. for i in range(numm)]

    mean = 1.
    var = 1.

    X=  [ [ random.gauss( mean,var) for f in range(n)  ]  for i in range(nump)] + \
        [ [ random.gauss(-mean,var) for f in range(n)  ]  for i in range(numm)] 

    primal_svm(m,n,X,y,CC)




qcqp_sdo_relaxation.py

Listing 16.34 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 *
import mosek.fusion.pythonic

# The relaxed SDP model
def miqcqp_sdo_relaxation(n,P,q):
    M = Model()

    M.setLogHandler(sys.stdout) 

    Z = M.variable("Z", Domain.inPSDCone(n+1))
    X = Z[:-1, :-1]
    x = Z[:-1, -1:]

    M.constraint( X.diag() >= x )
    M.constraint( Z[n, n] == 1 )

    M.objective( ObjectiveSense.Minimize, 
        Expr.dot( P, X ) + 2.0 * Expr.dot(x, q) )

    return M

# A direct integer model for minimizing |Ax-b|
def int_least_squares(n, A, b):
    M = Model()

    M.setLogHandler(sys.stdout) 

    x = M.variable("x", n, Domain.integral(Domain.unbounded()))
    t = M.variable("t", 1, Domain.unbounded())

    M.constraint( Expr.vstack(t, A @ x - b), Domain.inQCone() )
    M.objective( ObjectiveSense.Minimize, t )

    return M

# problem dimensions
n = 20
m = 2*n

# problem data
A = numpy.reshape(numpy.random.normal(0., 1.0, n*m), (m,n))
c = numpy.random.uniform(0., 1.0, n)
P = A.transpose().dot(A)
q = - P.dot(c)
b = A.dot(c)

# solve the problems
M = miqcqp_sdo_relaxation(n, P, q)
Mint = int_least_squares(n, A, b)
M.solve()
Mint.solve()

M.writeTask('M.ptf')
Mint.writeTask('Mint.ptf')

# rounded and optimal solution
xRound = numpy.rint(M.getVariable("Z").slice([0,n], [n,n+1]).level())
xOpt = numpy.rint(Mint.getVariable("x").level())

print(M.getSolverDoubleInfo("optimizerTime"), Mint.getSolverDoubleInfo("optimizerTime"))
print(numpy.linalg.norm(A.dot(xRound)-b), numpy.linalg.norm(A.dot(xOpt)-b)) 

reoptimization.py

Listing 16.35 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)
        print(x3.getShape())
        # Add to the exising constraint
        con.update(Expr.mul(x3, [4, 0, 1]), x3)
        # Change the objective to include x3
        M.objective(ObjectiveSense.Maximize, Expr.dot(c+[1.0], xNew))
        M.solve()
        print("x = {}".format(xNew.level()))

        ############## Add a new constraint ################
        con2 = M.constraint(Expr.dot(xNew, [1, 2, 1, 1]), Domain.lessThan(30000.0))
        M.solve()
        print("x = {}".format(xNew.level()))

        ############## Change constraint bounds ################
        cAll = Constraint.vstack(con, con2)
        # Change bounds by effectively updating fixed terms with the difference
        cAll.update([20000, 10000, 10000, 8000])
        M.solve()
        print("x = {}".format(xNew.level()))

response.py

Listing 16.36 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.37 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.38 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.39 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 mosek.fusion.pythonic
import sys
import numpy as np 

# Sample input data
n = 100
d = 4
k = 3
b = [9,10,11]
A = list(map(lambda x: x+np.transpose(x),
             [np.random.normal(0, 5, size=(d, d)) for _ in range(k*n)]))

# Create a model with n semidefinite variables od dimension d x d
with Model("sdo3") as M:
    X = M.variable(Domain.inPSDCone(d, n))

    # Pick indexes of diagonal entries for the objective
    alldiag = [[j,s,s] for j in range(n) for s in range(d)]
    M.objective(ObjectiveSense.Minimize, Expr.sum( X.pick(alldiag) ))

    # Each constraint is a sum of inner products
    # Each semidefinite variable is a slice of X
    for i in range(k):
        M.constraint(Expr.add([Expr.dot(A[i*n+j], X[j:j+1, :, :].reshape([d,d])) 
                               for j in range(n)]) >= b[i])

    # Solve
    M.setLogHandler(sys.stdout)            # Add logging
    M.writeTask("sdo3.ptf")                # Save problem in readable format
    M.solve()

    # Get results. Each variable is a slice of X
    print("Contributing variables:")
    for j in range(n):
        Xj = X.slice([j,0,0],[j+1,d,d]).level()
        if any(Xj[s]>1e-6 for s in range(d*d)):
            print("X{0}=\n{1}".format(j, np.reshape(Xj,(d,d))))

sospoly.py

Listing 16.40 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 *
import mosek.fusion.pythonic

def Hankel(n, i, a=1.0):
    '''Creates a Hankel matrix of dimension n+1, where 
       H_lk = a if l+k=i, and 0 otherwise.'''

    if i < n+1:
        return Matrix.sparse(n+1, n+1, range(i,-1,-1), range(i+1), (i+1)*[a])

    return Matrix.sparse(n+1, n+1, range(n, i-n-1, -1), range(i-n, n+1), (2*n+1-i)*[a])
    
def nn_inf(M, x):
    '''Models the cone of nonnegative polynomials on the real axis'''

    m = int(x.getSize() - 1)
    n = int(m/2) # degree of polynomial is 2*n

    assert(m == 2*n)    

    # Setup variables
    X = M.variable(Domain.inPSDCone(n+1))

    # x_i = Tr H(n, i) * X  i=0,...,m
    for i in range(m+1):
        M.constraint( x[i] == Expr.dot(Hankel(n,i), X) ) 

def nn_semiinf(M, x):
    '''Models the cone of nonnegative polynomials on the semi-infinite interval [0,inf)'''

    n = int(x.getSize()-1)
    n1, n2 = int(n/2), int((n-1)/2)
    
    # Setup variables
    X1 = M.variable(Domain.inPSDCone(n1+1))
    X2 = M.variable(Domain.inPSDCone(n2+1))

    # x_i = Tr H(n1, i) * X1 + Tr H(n2,i-1) * X2, i=0,...,n
    for i in range(n+1):       
        e1 = Expr.dot(Hankel(n1,i),  X1)
        e2 = Expr.dot(Hankel(n2,i-1),X2)        
        M.constraint( x[i] == e1 + e2 )
                 
def nn_finite(M, x, a, b):
    '''Models the cone of nonnegative polynomials on the finite interval [a,b]'''

    assert(a < b)
    m = int(x.getSize()-1)
    n = int(m/2)

    X1 = M.variable(Domain.inPSDCone(n+1))

    if (m == 2*n):
        X2 = M.variable(Domain.inPSDCone(n))

        # x_i = Tr H(n,i)*X1 + (a+b)*Tr H(n-1,i-1) * X2 - a*b*Tr H(n-1,i)*X2 - Tr H(n-1,i-2)*X2, i=0,...,m
        for i in range(m+1):        
            e1 = Expr.dot(Hankel(n  , i       ), X1) - Expr.dot(Hankel(n-1, i, a*b), X2)
            e2 = Expr.dot(Hankel(n-1, i-1, a+b), X2) - Expr.dot(Hankel(n-1, i-2   ), X2)
            M.constraint( x[i] == e1 + e2 )
    else:
        X2 = M.variable(Domain.inPSDCone(n+1))

        # x_i = Tr H(n,i-1)*X1 - a*Tr H(n,i)*X1 + b*Tr H(n,i)*X2 - Tr H(n,i-1)*X2, i=0,...,m
        for i in range(m+1):
            e1 = Expr.dot(Hankel(n, i-1 ), X1) - Expr.dot(Hankel(n, i, a), X1)
            e2 = Expr.dot(Hankel(n, i, b), X2) - Expr.dot(Hankel(n, i-1 ), X2)
            M.constraint( x[i] == e1 + e2 )

def diff(M, x):
    '''returns variables u representing the derivative of
    x(0) + x(1)*t + ... + x(n)*t^n,
    with u(0) = x(1), u(1) = 2*x(2), ..., u(n-1) = n*x(n).'''

    n = int(x.getSize()-1)
    u = M.variable(n, Domain.unbounded())

    mask = Matrix.dense(n,1,[float(i) for i in range(1,n+1)])

    M.constraint(u == Expr.mulElm(mask, x[1:n+1]))
    return u

def fitpoly(data, n):

    with Model('smooth poly') as M:

        m = len(data)
        A = [ [ data[j][0]**i for i in range(n+1) ] for j in range(m)]
        b = [ data[j][1] for j in range(m) ]
            
        x  = M.variable('x', n+1, Domain.unbounded())
        z  = M.variable('z', 1, Domain.unbounded())
        dx = diff(M, x)
                       
        M.constraint(A @ x == b)
                        
        # z - f'(t) >= 0, for all t \in [a, b]
        ub = M.variable(n, Domain.unbounded())

        dx0  = dx[0]
        dx1n = dx[1:n]

        M.constraint(ub == Expr.vstack(z - dx0, -dx1n))

        nn_finite(M, ub, data[0][0], data[-1][0])
            
        # f'(t) + z >= 0, for all t \in [a, b]
        lb = M.variable(n, Domain.unbounded())
        M.constraint(lb == Expr.vstack(z - dx0, dx1n))

        nn_finite(M, lb, data[0][0], data[-1][0])

        M.objective(ObjectiveSense.Minimize, z)
        M.solve()
        
        return x.level()

if __name__ == '__main__':

    data = [ [-1.0, 1.0], 
             [ 0.0, 0.0],
             [ 1.0, 1.0] ]
       
    x2 = fitpoly(data, 2)
    x4 = fitpoly(data, 4)
    x8 = fitpoly(data, 8)

    try:
        from pyx import *            
        from pyx.graph.axis import painter, tick

        text.set(mode="latex")

        p = painter.regular(basepathattrs=[deco.earrow.normal],
                            titlepos=0.95, 
                            titledist=-0.3,
                            titledirection=None,
                            outerticklength=graph.axis.painter.ticklength.normal,
                            )

        g = graph.graphxy(width=8, xaxisat=0, yaxisat=0,
                          x=graph.axis.linear(title="$t$", min=-2.9, max=2.9,
                                              painter=p),
                          y=graph.axis.linear(min=-0.05, max=1.9,
                                              manualticks=[tick.tick(0, None, None),
                                                           tick.tick(0.5, label=r"$\frac{1}{2}$", ticklevel=0),
                                                           tick.tick(1.5, label=r"$\frac{3}{2}$", ticklevel=0),
                                                           ],
                                              painter=p))
    
        def f(t,x): return t, sum([ x[i]*(t**i) for i in range(len(x)) ]) 

        g.plot(graph.data.paramfunction("t", -3, 3, "x, y = f(t,x)", points=500, context={"f": f, "x":x2}), 
               [graph.style.line([style.linewidth.Thin, style.linestyle.solid])])
        g.plot(graph.data.paramfunction("t", -3, 3, "x, y = f(t,x)", points=500, context={"f": f, "x":x4}),
               [graph.style.line([style.linewidth.Thin, style.linestyle.solid])])
        g.plot(graph.data.paramfunction("t", -3, 3, "x, y = f(t,x)", points=500, context={"f": f, "x":x8}),
               [graph.style.line([style.linewidth.Thin, style.linestyle.solid])])

        px, py = g.pos(1.3, f(1.3, x2)[1])
        g.text(px - 0.1, py, "$f_2(t)$", [text.halign.right, text.valign.top])
        
        px, py = g.pos(1.6, f(1.6, x4)[1])
        g.text(px + 0.1, py, "$f_4(t)$", [text.halign.left, text.valign.top])

        px, py = g.pos(1.31, f(1.31, x8)[1])
        g.text(px - 0.1, py, "$f_8(t)$", [text.halign.right, text.valign.top])

        g.writeEPSfile("sospoly")
        g.writePDFfile("sospoly")
        print("generated sospoly.eps")
    except ImportError:
        print("No PyX available")
        # No pyx available
        pass

sudoku.py

Listing 16.41 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 mosek.fusion.pythonic
import numpy as npy


def print_solution(m, x):
    n = m * m
    print("\n")
    for i in range(n):
        ss = ""
        for j in range(n):
            if j % m == 0:
                ss = ss + " |"

            for k in range(n):
                if x.index([i, j, k]).level() > 0.5:
                    ss = ss + " " + str(k + 1)
                    break

        print(ss + ' |')
        if (i + 1) % m == 0:
            print("\n")

def main():

    m = 3
    n = m * m

    hr_fixed = [[1, 5, 4],
                [2, 2, 5], [2, 3, 8], [2, 6, 3],
                [3, 2, 1], [3, 4, 2], [3, 5, 8], [3, 7, 9],
                [4, 2, 7], [4, 3, 3], [4, 4, 1], [4, 7, 8], [4, 8, 4],
                [6, 2, 4], [6, 3, 1], [6, 6, 9], [6, 7, 2], [6, 8, 7],
                [7, 3, 4], [7, 5, 6], [7, 6, 5], [7, 8, 8],
                [8, 4, 4], [8, 7, 1], [8, 8, 6],
                [9, 5, 9]
                ]

    fixed = [[f[0] - 1, f[1] - 1, f[2] - 1] for f in hr_fixed]

    with Model('SUDOKU') as M:
        x = M.variable("x",[n, n, n], Domain.binary())

        #each value only once per dimension
        for d in range(m):
            M.constraint("row_sum(%d)" % d, Expr.sum(x, d) == 1)

        #each number must appear only once in a block

        for k in range(n):
            for i in range(m):
                for j in range(m):
                    M.constraint("blocksum(%d,%d,k=%d)" % (i,j,k),
                                 Expr.sum(x[i*m : (i+1)*m, j*m : (j+1)*m, k : (k+1)]) == 1)

        M.constraint("fix", x.pick(fixed) == 1.0)

        M.setLogHandler(sys.stdout)

        M.solve()

        M.writeTask("sudoku.task")

        #print the solution, if any...
        if M.getPrimalSolutionStatus() in [SolutionStatus.Optimal]:
            print_solution(m, x)
        else:
            print("No solution found!")

if __name__ == '__main__':
    main()

total_variation.py

Listing 16.42 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 mosek.fusion.pythonic
import numpy as np

def total_var(n,m):
    M = Model('TV')

    u = M.variable("u", [n+1,m+1], Domain.inRange(0.,1.0) )
    t = M.variable("t", [n,m], Domain.unbounded() )

    # In this example we define sigma and the input image f as parameters
    # to demonstrate how to solve the same model with many data variants.
    # Of course they could simply be passed as ordinary arrays if that is not needed.
    sigma = M.parameter("sigma")
    f = M.parameter("f", [n,m])

    ucore = u[:-1, :-1]

    deltax = u[1:, :-1] - ucore
    deltay = u[:-1, 1:] - ucore

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

    M.constraint( Expr.vstack(sigma, Expr.flatten( f - ucore )),
                  Domain.inQCone() )

    M.objective( ObjectiveSense.Minimize, Expr.sum(t) )

    return M

#Display
def show(n,m,grid):
    try:
        import matplotlib
        import matplotlib.pyplot as plt
        import matplotlib.cm as cm
        plt.imshow(grid, extent=(0,m,0,n),
                   interpolation='nearest', cmap=cm.jet,
                   vmin=0, vmax=1)
        plt.show()
    except:
        print (grid)

if __name__ == '__main__':
    np.random.seed(0)
    
    n, m = 100, 200

    # Create a parametrized model with given shape
    M = total_var(n, m)
    sigma = M.getParameter("sigma")
    f     = M.getParameter("f")
    ucore = M.getVariable("u").slice([0,0], [n,m])

    # Example: Linear signal with Gaussian noise    
    signal = np.reshape([[1.0*(i+j)/(n+m) for i in range(m)] for j in range(n)], (n,m))
    noise  = np.random.normal(0., 0.08, (n,m))
    fVal   = signal + noise

    # Uncomment to get graphics:
    # show(n, m, signal)
    # show(n, m, fVal)
   
    # Set value for f
    f.setValue(fVal)

    for sigmaVal in [0.0004, 0.0005, 0.0006]:
        # Set new value for sigma and solve
        sigma.setValue(sigmaVal*n*m)

        M.solve()

        sol = np.reshape(ucore.level(), (n,m))
        # Now use the solution
        # ...

        # Uncomment to get graphics:
        # show(n, m, np.reshape(ucore.level(), (n,m)))

        print("rel_sigma = {sigmaVal}  total_var = {var}".format(sigmaVal=sigmaVal,
                                                                 var=M.primalObjValue() ))


    M.dispose()

tsp.py

Listing 16.43 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 mosek.fusion.pythonic
import sys

def tsp(n, A, C, remove_selfloops, remove_2_hop_loops):
    with Model() as M:

        M.setLogHandler(sys.stdout)
        x = M.variable([n,n], Domain.binary())

        M.constraint(Expr.sum(x,0) == 1.0)
        M.constraint(Expr.sum(x,1) == 1.0)
        M.constraint(x <= A)

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

        if remove_2_hop_loops:
            M.constraint(x + x.T <= 1.0)

        if remove_selfloops:
            M.constraint(x.diag() == 0)

        it = 1
        M.writeTask("tsp-0-%s-%s.ptf" % ('t' if remove_selfloops else 'f', 't' if remove_2_hop_loops else 'f'))

        while True:
            print("\n\n--------------------\nIteration",it)
            M.solve()

            print('\nsolution cost:', M.primalObjValue())
            print('\nsolution:')

            cycles = []

            for i in range(n):
                xi = x[i:i+1, :]
                print(xi.level())

                for j in range(n):
                    if xi.level()[j] <= 0.5 : continue

                    found = False
                    for c in cycles:
                        if len( [ a for a in c if i in a or j in a ] )> 0:
                            c.append( [i,j] )
                            found = True
                            break

                    if not found:
                        cycles.append([ [ i,j ]])

            print('\ncycles:')
            print([c for c in cycles])

            if len(cycles)==1:
                break;

            for c in cycles:
                M.constraint(Expr.sum(x.pick(c)) <= len(c) - 1)
            it = it +1

        return x.level(), c

    return [],[]

def main():
    A_i = [0,1,2,3,1,0,2,0]
    A_j = [1,2,3,0,0,2,1,3]
    C_v = [1.,1.,1.,1.,0.1,0.1,0.1,0.1]
    n = max(max(A_i),max(A_j))+1
    costs = Matrix.sparse(n,n,A_i,A_j,C_v)
    x,c = tsp(n, Matrix.sparse(n,n,A_i,A_j,1.), costs , True, True)
    x,c = tsp(n, Matrix.sparse(n,n,A_i,A_j,1.), costs , True, False)
if __name__ == '__main__':
    main()