17 List of examples

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

Table 17.1 List of distributed examples

File

Description

acc1.py

A simple problem with one affine conic constraint (ACC)

acc2.py

A simple problem with two affine conic constraints (ACC)

blas_lapack.py

Demonstrates the MOSEK interface to BLAS/LAPACK linear algebra routines

callback.py

An example of data/progress callback

ceo1.py

A simple conic exponential problem

concurrent1.py

Implementation of a concurrent optimizer for linear and mixed-integer problems

cqo1.py

A simple conic quadratic problem

djc1.py

A simple problem with disjunctive constraints (DJC)

feasrepairex1.py

A simple example of how to repair an infeasible problem

gp1.py

A simple geometric program (GP) in conic form

helloworld.py

A Hello World example

lo1.py

A simple linear problem

lo2.py

A simple linear problem

logistic.py

Implements logistic regression and simple log-sum-exp (CEO)

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

opt_server_async.py

Uses MOSEK OptServer to solve an optimization problem asynchronously

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

qcqo1.py

A simple quadratically constrained quadratic problem

qo1.py

A simple quadratic 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

sdo_lmi.py

A simple semidefinite problem with an LMI using the SVEC domain.

sensitivity.py

Sensitivity analysis performed on a small linear problem

simple.py

A simple I/O example: read problem from a file, solve and write solutions

solutionquality.py

Demonstrates how to examine the quality of a solution

solvebasis.py

Demonstrates solving a linear system with the basis matrix

solvelinear.py

Demonstrates solving a general linear system

sparsecholesky.py

Shows how to find a Cholesky factorization of a sparse matrix

Additional examples can be found on the MOSEK website and in other MOSEK publications.

acc1.py

Listing 17.1 acc1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      acc1.py
#
#  Purpose :   Tutorial example for affine conic constraints.
#              Models the problem:
#
#              maximize c^T x
#              subject to  sum(x) = 1
#                          gamma >= |Gx+h|_2
##
import sys, mosek

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()

# Define problem data
n, k = 3, 2
# Only a symbolic constant
inf = 0.0

# Make a MOSEK environment
with mosek.Env() as env:
    # Attach a printer to the environment
    env.set_Stream(mosek.streamtype.log, streamprinter)

    # Create a task
    with env.Task(0, 0) as task:
        # Attach a printer to the task
        task.set_Stream(mosek.streamtype.log, streamprinter)

        # Create n free variables
        task.appendvars(n)
        task.putvarboundsliceconst(0, n, mosek.boundkey.fr, -inf, inf)

        # Set up the objective
        c = [2, 3, -1]
        task.putobjsense(mosek.objsense.maximize)
        task.putclist(range(n), c)

        # One linear constraint - sum(x) = 1
        task.appendcons(1)
        task.putarow(0, range(n), [1]*n)
        task.putconbound(0, mosek.boundkey.fx, 1.0, 1.0)

        # Append empty AFE rows for affine expression storage
        task.appendafes(k + 1)

        # G matrix in sparse form
        Gsubi = [0, 0, 1, 1]
        Gsubj = [0, 1, 0, 2]
        Gval  = [1.5, 0.1, 0.3, 2.1]
        # Other data
        h     = [0, 0.1]
        gamma = 0.03

        # Construct F matrix in sparse form
        Fsubi = [i + 1 for i in Gsubi]   # G will be placed from row number 1 in F
        Fsubj = Gsubj
        Fval  = Gval

        # Fill in F storage
        task.putafefentrylist(Fsubi, Fsubj, Fval)

        # Fill in g storage
        task.putafeg(0, gamma)
        task.putafegslice(1, k+1, h)

        # Define a conic quadratic domain
        quadDom = task.appendquadraticconedomain(k + 1)

        # Create the ACC
        task.appendacc(quadDom,    # Domain index
                       range(k+1), # Indices of AFE rows [0,...,k]
                       None)       # Ignored

        task.writedata("acc1.ptf")
        # Solve and retrieve solution
        task.optimize()
        xx = task.getxx(mosek.soltype.itr)
        if task.getsolsta(mosek.soltype.itr) == mosek.solsta.optimal:
            print("Solution: {xx}".format(xx=list(xx)))

        # Demonstrate retrieving activity of ACC
        activity = task.evaluateacc(mosek.soltype.itr, 
                                    0)     # ACC index          
        print("Activity of ACC:: {activity}".format(activity=list(activity)))

        # Demonstrate retrieving the dual of ACC
        doty = task.getaccdoty(mosek.soltype.itr, 
                               0)          # ACC index
        print("Dual of ACC:: {doty}".format(doty=list(doty)))

acc2.py

Listing 17.2 acc2.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      acc2.py
#
#  Purpose :   Tutorial example for affine conic constraints.
#              Models the problem:
#
#              maximize c^T x
#              subject to  sum(x) = 1
#                          gamma >= |Gx+h|_2
#
#              This version inputs the linear constraint as an affine conic constraint.
##
import sys, mosek

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()

# Define problem data
n, k = 3, 2
# Only a symbolic constant
inf = 0.0

# Make a MOSEK environment
with mosek.Env() as env:
    # Attach a printer to the environment
    env.set_Stream(mosek.streamtype.log, streamprinter)

    # Create a task
    with env.Task(0, 0) as task:
        # Attach a printer to the task
        task.set_Stream(mosek.streamtype.log, streamprinter)

        # Create n free variables
        task.appendvars(n)
        task.putvarboundsliceconst(0, n, mosek.boundkey.fr, -inf, inf)

        # Set up the objective
        c = [2, 3, -1]
        task.putobjsense(mosek.objsense.maximize)
        task.putclist(range(n), c)

        # Set AFE rows representing the linear constraint
        task.appendafes(1)
        task.putafefrow(0, range(n), [1]*n)
        task.putafeg(0, -1.0)

        # Set AFE rows representing the quadratic constraint
        task.appendafes(k + 1)
        task.putafefrow(2,          # afeidx, row number
                        [0, 1],     # varidx, column numbers
                        [1.5, 0.1]) # values
        task.putafefrow(3,          # afeidx, row number
                        [0, 2],     # varidx, column numbers
                        [0.3, 2.1]) # values

        h     = [0, 0.1]
        gamma = 0.03
        task.putafeg(1, gamma)
        task.putafegslice(2, k+2, h)

        # Define domains
        zeroDom = task.appendrzerodomain(1)
        quadDom = task.appendquadraticconedomain(k + 1)

        # Append affine conic constraints
        task.appendacc(zeroDom,    # Domain index
                       [0],        # Indices of AFE rows
                       None)       # Ignored
        task.appendacc(quadDom,    # Domain index
                       [1,2,3],    # Indices of AFE rows
                       None)       # Ignored

        # Solve and retrieve solution
        task.optimize()
        xx = task.getxx(mosek.soltype.itr)
        if task.getsolsta(mosek.soltype.itr) == mosek.solsta.optimal:
            print("Solution: {xx}".format(xx=list(xx)))

        # Demonstrate retrieving activity of ACC
        activity = task.evaluateacc(mosek.soltype.itr, 
                                    1)     # ACC index
        print("Activity of quadratic ACC:: {activity}".format(activity=list(activity)))

        # Demonstrate retrieving the dual of ACC
        doty = task.getaccdoty(mosek.soltype.itr, 
                               1)          # ACC index
        print("Dual of quadratic ACC:: {doty}".format(doty=list(doty)))

blas_lapack.py

Listing 17.3 blas_lapack.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      blas_lapack.py
#
#  Purpose :   To demonstrate how to call BLAS/LAPACK routines
#              for whose MOSEK provides simplified interfaces.
##
import mosek

def print_matrix(x, r, c):
    for i in range(r):
        print([x[j * r + i] for j in range(c)])

with mosek.Env() as env:

    n = 3
    m = 2
    k = 3

    alpha = 2.0
    beta = 0.5

    x = [1.0, 1.0, 1.0]
    y = [1.0, 2.0, 3.0]
    z = [1.0, 1.0]
    v = [0.0, 0.0]
    #A has m=2 rows and k=3 cols
    A = [1.0, 1.0, 2.0, 2.0, 3., 3.]
    #B has k=3 rows and n=3 cols
    B = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    C = [0.0 for i in range(n * m)]
    D = [1.0, 1.0, 1.0, 1.0]
    Q = [1.0, 0.0, 0.0, 2.0]

# BLAS routines

    xy = env.dot(n, x, y)
    print("dot results= %f\n" % xy)

    env.axpy(n, alpha, x, y)
    print("\naxpy results is ")
    print_matrix(y, 1, len(y))

    env.gemv(mosek.transpose.no, m, n, alpha, A, x, beta, z)
    print("\ngemv results is ")
    print_matrix(z, 1, len(z))

    env.gemm(mosek.transpose.no, mosek.transpose.no,
             m, n, k, alpha, A, B, beta, C)
    print("\ngemm results is ")
    print_matrix(C, m, n)

    env.syrk(mosek.uplo.lo, mosek.transpose.no, m, k, alpha, A, beta, D)
    print("\nsyrk results is")
    print_matrix(D, m, m)

# LAPACK routines

    env.potrf(mosek.uplo.lo, m, Q)
    print("\npotrf results is ")
    print_matrix(Q, m, m)


    env.syeig(mosek.uplo.lo, m, Q, v)
    print("\nsyeig results is")
    print_matrix(v, 1, m)


    env.syevd(mosek.uplo.lo, m, Q, v)
    print("\nsyevd results is")
    print('v: ')
    print_matrix(v, 1, m)
    print('Q: ')
    print_matrix(Q, m, m)

    print("Exiting...")

callback.py

Listing 17.4 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.
#
#              Use this script as follows:
#
#               callback.py psim  25fv47.mps
#               callback.py dsim  25fv47.mps
#               callback.py intpnt 25fv47.mps
#
#              The first argument tells which optimizer to use
#              i.e. psim is primal simplex, dsim is dual simplex
#              and intpnt is interior-point.
##
import sys
import numpy

import mosek
from mosek import *

def makeUserCallback(maxtime, task):
    pass

    def userCallback(caller,
                     douinf,
                     intinf,
                     lintinf):
        opttime = 0.0

        if caller == callbackcode.begin_intpnt:
            print("Starting interior-point optimizer")
        elif caller == callbackcode.intpnt:
            itrn = intinf[iinfitem.intpnt_iter]
            pobj = douinf[dinfitem.intpnt_primal_obj]
            dobj = douinf[dinfitem.intpnt_dual_obj]
            stime = douinf[dinfitem.intpnt_time]
            opttime = douinf[dinfitem.optimizer_time]

            print("Iterations: %-3d" % itrn)
            print("  Elapsed time: %6.2f(%.2f) " % (opttime, stime))
            print("  Primal obj.: %-18.6e  Dual obj.: %-18.6e" % (pobj, dobj))
        elif caller == callbackcode.end_intpnt:
            print("Interior-point optimizer finished.")
        elif caller == callbackcode.begin_primal_simplex:
            print("Primal simplex optimizer started.")
        elif caller == callbackcode.update_primal_simplex:
            itrn = intinf[iinfitem.sim_primal_iter]
            pobj = douinf[dinfitem.sim_obj]
            stime = douinf[dinfitem.sim_time]
            opttime = douinf[dinfitem.optimizer_time]

            print("Iterations: %-3d" % itrn)
            print("  Elapsed time: %6.2f(%.2f)" % (opttime, stime))
            print("  Obj.: %-18.6e" % pobj)
        elif caller == callbackcode.end_primal_simplex:
            print("Primal simplex optimizer finished.")
        elif caller == callbackcode.begin_dual_simplex:
            print("Dual simplex optimizer started.")
        elif caller == callbackcode.update_dual_simplex:
            itrn = intinf[iinfitem.sim_dual_iter]
            pobj = douinf[dinfitem.sim_obj]
            stime = douinf[dinfitem.sim_time]
            opttime = douinf[dinfitem.optimizer_time]
            print("Iterations: %-3d" % itrn)
            print("  Elapsed time: %6.2f(%.2f)" % (opttime, stime))
            print("  Obj.: %-18.6e" % pobj)
        elif caller == callbackcode.end_dual_simplex:
            print("Dual simplex optimizer finished.")
        elif caller == callbackcode.new_int_mio:
            print("New integer solution has been located.")
            xx = task.getxx(soltype.itg)
            print(xx)
            print("Obj.: %f" % douinf[dinfitem.mio_obj_int])
        else:
            pass

        if opttime >= maxtime:
            # mosek is spending too much time. Terminate it.
            print("Terminating.")
            return 1

        return 0
    return userCallback


def main(args):
    # To run a continuous problem example
    filename = "../data/25fv47.mps"
    slvr = "intpnt"

    # To run a mixed-integer example
    #filename = "../data/milo1.lp"
    #slvr     = ""

    if len(args) < 3:
        print("Usage: callback ( psim | dsim | intpnt ) filename")

    if len(args) > 1:
        slvr = args[1]
    if len(args) > 2:
        filename = args[2]

    with mosek.Env() as env:
        with mosek.Task(env) as task:
            task.readdata(filename)

            if slvr == 'psim':
                task.putintparam(iparam.optimizer,
                                 optimizertype.primal_simplex)
            elif slvr == "dsim":
                task.putintparam(iparam.optimizer, optimizertype.dual_simplex)
            elif slvr == "intpnt":
                task.putintparam(iparam.optimizer, optimizertype.intpnt)

            usercallback = makeUserCallback(maxtime=0.05, task=task)
            task.set_InfoCallback(usercallback)

            task.optimize()


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

ceo1.py

Listing 17.5 ceo1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      ceo1.py
#
#  Purpose :   Demonstrates how to solve small conic exponential
#              optimization problem using the MOSEK Python API.
##
import sys
import mosek

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()

def main():

    # Only a symbolic constant
    inf = 0.0

    # Create a task
    with mosek.Task() as task:
        # Attach a printer to the task
        task.set_Stream(mosek.streamtype.log, streamprinter)

        c = [1.0, 1.0, 0.0]
        a = [1.0, 1.0, 1.0]
        numvar, numcon = 3, 1

        # Append 'numcon' empty constraints.
        # The constraints will initially have no bounds.
        task.appendcons(numcon)

        # Append 'numvar' variables.
        # The variables will initially be fixed at zero (x=0).
        task.appendvars(numvar)

        # Set up the linear part of the problem
        task.putcslice(0, numvar, c)
        task.putarow(0, [0, 1, 2], a)
        task.putvarboundslice(0, numvar, [mosek.boundkey.fr] * numvar, [inf] * numvar, [inf] * numvar)
        task.putconbound(0, mosek.boundkey.fx, 1.0, 1.0)

        # Add a conic constraint
        # Create a 3x3 identity matrix F
        task.appendafes(3)
        task.putafefentrylist([0, 1, 2],         # Rows 
                              [0, 1, 2],         # Columns 
                              [1.0] * 3)

        # Exponential cone (x(0),x(1),x(2)) \in EXP 
        expdomain  = task.appendprimalexpconedomain()
        task.appendacc(expdomain,               # Domain
                       [0, 1, 2],               # Rows from F 
                       None)                    # Unused 

        # Input the objective sense (minimize/maximize)
        task.putobjsense(mosek.objsense.minimize)

        # Optimize the task
        task.optimize()
        # Print a summary containing information
        # about the solution for debugging purposes
        task.solutionsummary(mosek.streamtype.msg)
        prosta = task.getprosta(mosek.soltype.itr)
        solsta = task.getsolsta(mosek.soltype.itr)

        # Output a solution
        xx = task.getxx(mosek.soltype.itr)

        if solsta == mosek.solsta.optimal:
            print("Optimal solution: %s" % xx)
        elif solsta == mosek.solsta.dual_infeas_cer:
            print("Primal or dual infeasibility.\n")
        elif solsta == mosek.solsta.prim_infeas_cer:
            print("Primal or dual infeasibility.\n")
        elif mosek.solsta.unknown:
            print("Unknown solution status")
        else:
            print("Other solution status")


# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.code))
    if msg is not None:
        print("\t%s" % e.msg)
        sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

concurrent1.py

Listing 17.6 concurrent1.py Click here to download.
#
#  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File:      concurrent1.py
#
#  Purpose: Demonstrates a simple implementation of a concurrent optimizer.
#
#           The concurrent optimizer starts a few parallel optimizations
#           of the same problem using different algorithms, and reports
#           a solution when the first optimizer is ready.
#
#           This example also demonstrates how to define a simple callback handler
#           that stops the optimizer when requested.
#
import mosek, sys
from threading import Thread

# Defines a Mosek callback function whose only function
# is to indicate if the optimizer should be stopped.
stop = False
firstStop = -1
def cbFun(code):
  return 1 if stop else 0

# Stream handling method
def streamwriter(s):
  sys.stdout.write('{0}'.format(s))

# firstOK, res, trm = optimize(tasks)
#
# Takes a list of tasks and optimizes then in parallel. The
# response code and termination code from each optimization is
# returned in ``res`` and ``trm``.
#
# When one task completes with rescode == ok, others are terminated.
#
# firstOK is the index of the first task that returned
# with rescode == ok. Whether or not this task contains the
# most valuable answer, is for the caller to decide. If none
# completed without error returns -1.
def runTask(num, task, res, trm):
  global stop
  global firstStop
  try:
    trm[num] = task.optimize();
    res[num] = mosek.rescode.ok
  except mosek.MosekException as e:
    trm[num] = mosek.rescode.err_unknown
    res[num] = e.errno
  finally:
    # If this finished with success, inform other tasks to interrupt
    if res[num] == mosek.rescode.ok:
      if not stop:
        firstStop = num
      stop = True

def optimize(tasks):
  n = len(tasks)
  res = [ mosek.rescode.err_unknown ] * n
  trm = [ mosek.rescode.err_unknown ] * n  

  # Set a callback function 
  for t in tasks:
    t.set_Progress(cbFun)
  
  # Start parallel optimizations, one per task
  jobs = [ Thread(target=runTask, args=(i, tasks[i], res, trm)) for i in range(n) ]
  for j in jobs:
    j.start()
  for j in jobs:
    j.join()

  # For debugging, print res and trm codes for all optimizers
  for i in range(n):
    print("Optimizer  {0}   res {1}   trm {2}".format(i, res[i], trm[i]))

  return firstStop, res, trm

#
# idx, winTask, winTrm, winRes = optimizeconcurrent(task, optimizers)
#
# Given a continuous task, set up jobs to optimize it 
# with a list of different solvers.
#
# Returns an index, corresponding to the optimization
# task that is returned as winTask. This is the task
# with the best possible status of those that finished.
# If none task is considered successful returns -1.
def optimizeconcurrent(task, optimizers):
  n = len(optimizers)
  tasks = [ mosek.Task(task) for _ in range(n) ]

  # Choose various optimizers for cloned tasks
  for i in range(n):
    tasks[i].putintparam(mosek.iparam.optimizer, optimizers[i])

  # Solve tasks in parallel
  firstOK, res, trm = optimize(tasks)

  if firstOK >= 0:
    return firstOK, tasks[firstOK], trm[firstOK], res[firstOK]
  else:
    return -1, None, None, None

#
# idx, winTask, winTrm, winRes = optimizeconcurrent(task, optimizers)
#
# Given a mixed-integer task, set up jobs to optimize it 
# with different values of seed. That will lead to
# different execution paths of the optimizer.
#
# Returns an index, corresponding to the optimization
# task that is returned as winTask. This is the task
# with the best value of the objective function.
# If none task is considered successful returns -1.
#
# Typically, the input task would contain a time limit. The two
# major scenarios are:
# 1. Some clone ends before time limit - then it has optimum.
# 2. All clones reach time limit - pick the one with best objective.
def optimizeconcurrentMIO(task, seeds):
  n = len(seeds)
  tasks = [ mosek.Task(task) for _ in range(n) ]  
  
  # Choose various seeds for cloned tasks
  for i in range(n):
    tasks[i].putintparam(mosek.iparam.mio_seed, seeds[i])

  # Solve tasks in parallel
  firstOK, res, trm = optimize(tasks)

  if firstOK >= 0:
    # Pick the task that ended with res = ok
    # and contains an integer solution with best objective value
    sense = task.getobjsense();
    bestObj = 1.0e+10 if sense == mosek.objsense.minimize else -1.0e+10
    bestPos = -1

    for i in range(n):
      print("{0}   {1}".format(i,tasks[i].getprimalobj(mosek.soltype.itg)))

    for i in range(n):
      if ((res[i] == mosek.rescode.ok) and
          (tasks[i].getsolsta(mosek.soltype.itg) == mosek.solsta.prim_feas or
           tasks[i].getsolsta(mosek.soltype.itg) == mosek.solsta.integer_optimal) and
         ((tasks[i].getprimalobj(mosek.soltype.itg) < bestObj) 
           if (sense == mosek.objsense.minimize) else 
          (tasks[i].getprimalobj(mosek.soltype.itg) > bestObj))):
        bestObj = tasks[i].getprimalobj(mosek.soltype.itg)
        bestPos = i

    if bestPos >= 0:
      return bestPos, tasks[bestPos], trm[bestPos], res[bestPos]

  return -1, None, None, None

 
# This is an example of how one can use the methods
#       optimizeconcurrent
#       optimizeconcurrentMIO
#
#   argv[0] : name of file with input problem
#   argv[1]: (optional) time limit
def main(argv):
  with mosek.Env() as env:
    with mosek.Task(env, 0, 0) as task:
      if len(argv) >= 2:
        task.readdata(argv[1])
      else:
        task.readdata("../data/25fv47.mps")

      # Optional time limit
      if len(argv) >= 3:
        task.putdouparam(mosek.dparam.optimizer_max_time, float(argv[2]))

      if (task.getnumintvar() == 0):
        # If the problem is continuous
        # optimize it with three continuous optimizers.
        # (Simplex will fail for non-linear problems)
        optimizers = [
          mosek.optimizertype.conic,
          mosek.optimizertype.dual_simplex,
          mosek.optimizertype.primal_simplex
        ]

        idx, t, trm, res = optimizeconcurrent(task, optimizers)
      else:
        # Mixed-integer problem.
        # Try various seeds.
        seeds = [ 42, 13, 71749373 ]

        idx, t, trm, res = optimizeconcurrentMIO(task, seeds)

      # Check results and print the best answer
      if idx >= 0:
        print("Result from optimizer with index {0}:  res {1}  trm {2}".format(idx, res, trm))
        t.set_Stream(mosek.streamtype.log, streamwriter)
        t.optimizersummary(mosek.streamtype.log)
        t.solutionsummary(mosek.streamtype.log);
      else:
        print("All optimizers failed.")


if __name__ == "__main__":
    main(sys.argv)

cqo1.py

Listing 17.7 cqo1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      cqo1.py
#
#  Purpose :   Demonstrates how to solve small conic quadratic
#              optimization problem using the MOSEK Python API.
##
import sys
import mosek

# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Create a task
    with mosek.Task() as task:
        # Attach a printer to the task
        task.set_Stream(mosek.streamtype.log, streamprinter)

        bkc = [mosek.boundkey.fx]
        blc = [1.0]
        buc = [1.0]

        c = [0.0, 0.0, 0.0,
             1.0, 1.0, 1.0]
        bkx = [mosek.boundkey.lo, mosek.boundkey.lo, mosek.boundkey.lo,
               mosek.boundkey.fr, mosek.boundkey.fr, mosek.boundkey.fr]
        blx = [0.0, 0.0, 0.0,
               -inf, -inf, -inf]
        bux = [inf, inf, inf,
               inf, inf, inf]

        asub = [[0], [0], [0]]
        aval = [[1.0], [1.0], [2.0]]

        numvar = len(bkx)
        numcon = len(bkc)
        NUMANZ = 4

        # Append 'numcon' empty constraints.
        # The constraints will initially have no bounds.
        task.appendcons(numcon)

        #Append 'numvar' variables.
        # The variables will initially be fixed at zero (x=0).
        task.appendvars(numvar)

        for j in range(numvar):
            # Set the linear term c_j in the objective.
            task.putcj(j, c[j])
            # Set the bounds on variable j
            # blx[j] <= x_j <= bux[j]
            task.putvarbound(j, bkx[j], blx[j], bux[j])

        for j in range(len(aval)):
            # Input column j of A
            task.putacol(j,                  # Variable (column) index.
                         # Row index of non-zeros in column j.
                         asub[j],
                         aval[j])            # Non-zero Values of column j.

        for i in range(numcon):
            task.putconbound(i, bkc[i], blc[i], buc[i])

        # Input the affine conic constraints
        # Create a matrix F such that F * x = [x(3),x(0),x(1),x(4),x(5),x(2)] 
        task.appendafes(6)
        task.putafefentrylist(range(6),                      # Rows
                              [3, 0, 1, 4, 5, 2],            # Columns 
                              [1.0] * 6)

        # Quadratic cone (x(3),x(0),x(1)) \in QUAD_3 
        quadcone  = task.appendquadraticconedomain(3)
        task.appendacc(quadcone,                # Domain 
                       [0, 1, 2],               # Rows from F
                       None)                    # Unused

        # Rotated quadratic cone (x(4),x(5),x(2)) \in RQUAD_3 
        rquadcone = task.appendrquadraticconedomain(3);
        task.appendacc(rquadcone,               # Domain 
                       [3, 4, 5],               # Rows from F 
                       None)                    # Unused 

        # Input the objective sense (minimize/maximize)
        task.putobjsense(mosek.objsense.minimize)

        # Optimize the task
        task.optimize()
        task.writedata("cqo1.ptf")
        # Print a summary containing information
        # about the solution for debugging purposes
        task.solutionsummary(mosek.streamtype.msg)
        prosta = task.getprosta(mosek.soltype.itr)
        solsta = task.getsolsta(mosek.soltype.itr)

        # Output a solution
        xx = task.getxx(mosek.soltype.itr)

        if solsta == mosek.solsta.optimal:
            print("Optimal solution: %s" % xx)
        elif solsta == mosek.solsta.dual_infeas_cer:
            print("Primal or dual infeasibility.\n")
        elif solsta == mosek.solsta.prim_infeas_cer:
            print("Primal or dual infeasibility.\n")
        elif mosek.solsta.unknown:
            print("Unknown solution status")
        else:
            print("Other solution status")


# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    print("\t%s" % e.msg)
    sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

djc1.py

Listing 17.8 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
##
import sys
from mosek import *

# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0

def main():
    # Make mosek environment
    with Env() as env:
        # Create a task object
        with env.Task(0, 0) as task:
            # Append free variables
            numvar = 4
            task.appendvars(numvar)
            task.putvarboundsliceconst(0, numvar, boundkey.fr, -inf, inf)

            # The linear part: the linear constraint
            task.appendcons(1)
            task.putarow(0, range(numvar), [1] * numvar)
            task.putconbound(0, boundkey.lo, -10.0, -10.0)

            # The linear part: objective
            task.putobjsense(objsense.minimize)
            task.putclist(range(numvar), [2, 1, 3, 1])

            # Fill in the affine expression storage F, g
            numafe = 10
            task.appendafes(numafe)

            fafeidx = [0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9]
            fvaridx = [0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]
            fval    = [1.0, -2.0, 1.0, -3.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
            g       = [1.0, 2.0, 0.0, 0.0, 0.0, 0.0, -2.5, -2.5, -2.5, -2.5]

            task.putafefentrylist(fafeidx, fvaridx, fval)
            task.putafegslice(0, numafe, g)

            # Create domains
            zero1   = task.appendrzerodomain(1)
            zero2   = task.appendrzerodomain(2)
            rminus1 = task.appendrminusdomain(1)

            # Append disjunctive constraints
            numdjc = 2
            task.appenddjcs(numdjc)

            # First disjunctive constraint
            task.putdjc(0,                                        # DJC index
                        [rminus1, zero2, rminus1, zero2],         # Domains     (domidxlist)
                        [0, 4, 5, 1, 2, 3],                       # AFE indices (afeidxlist)
                        None,                                     # Unused
                        [2, 2] )                                  # Term sizes  (termsizelist)

            # Second disjunctive constraint
            task.putdjc(1,                                        # DJC index
                        [zero1, zero1, zero1, zero1],             # Domains     (domidxlist)
                        [6, 7, 8, 9],                             # AFE indices (afeidxlist)
                        None,                                     # Unused
                        [1, 1, 1, 1] )                            # Term sizes  (termidxlist)

            # Useful for debugging
            task.writedata("djc.ptf")                          # Write file in human-readable format
            task.set_Stream(streamtype.log, sys.stdout.write)  # Attach a log stream printer to the task

            # Solve the problem
            task.optimize()

            # Print a summary containing information
            # about the solution for debugging purposes
            task.solutionsummary(streamtype.msg)

            # Get status information about the solution
            sta = task.getsolsta(soltype.itg)

            if (sta == solsta.integer_optimal):
                xx = task.getxx(soltype.itg)
                
                print("Optimal solution: ")
                for i in range(numvar):
                    print("x[" + str(i) + "]=" + str(xx[i]))
            else:
                print("Another solution status")

# call the main function
try:
    main()
except Error as e:
    print("ERROR: %s" % str(e.errno))
    if e.msg is not None:
        print("\t%s" % e.msg)
        sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

feasrepairex1.py

Listing 17.9 feasrepairex1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      feasrepairex1.py
#
#  Purpose :   To demonstrate how to use the MSK_relaxprimal function to
#              locate the cause of an infeasibility.
#
#  Syntax :    On command line
#              python feasrepairex1.py feasrepair.lp
#              feasrepair.lp is located in mosek\<version>\tools\examples.
##
import sys
import mosek

# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main(inputfile):
    # Make a MOSEK environment
    with mosek.Env() as env:
        with env.Task(0, 0) as task:
            # Attach a printer to the task
            task.set_Stream(mosek.streamtype.log, streamprinter)

            # Read data
            task.readdata(inputfile)

            task.putintparam(mosek.iparam.log_feas_repair, 3)

            task.primalrepair(None, None, None, None)

            sum_viol = task.getdouinf(mosek.dinfitem.primal_repair_penalty_obj)
            print("Minimized sum of violations = %e" % sum_viol)

            task.optimize()

            task.solutionsummary(mosek.streamtype.msg)

# call the main function
try:
    filename = "../data/feasrepair.lp"
    if len(sys.argv) > 1:
        filename = sys.argv[1]
    main(filename)
except Exception as e:
    print(e)
    raise

gp1.py

Listing 17.10 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
from mosek import *
import sys

# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()

# maximize     h*w*d
# subjecto to  2*(h*w + h*d) <= Awall
#              w*d <= Afloor
#              alpha <= h/w <= beta
#              gamma <= d/w <= delta
#
# Variable substitutions:  h = exp(x), w = exp(y), d = exp(z).
#
# maximize     x+y+z
# subject      log( exp(x+y+log(2/Awall)) + exp(x+z+log(2/Awall)) ) <= 0
#                              y+z <= log(Afloor)
#              log( alpha ) <= x-y <= log( beta )
#              log( gamma ) <= z-y <= log( delta )
def max_volume_box(Aw, Af, alpha, beta, gamma, delta):
    # Basic dimensions of our problem
    numvar    = 3       # Variables in original problem
    x, y, z   = 0, 1, 2 # Indices of variables
    numcon    = 3       # Linear constraints in original problem

    # Linear part of the problem
    cval  = [1, 1, 1]
    asubi = [0, 0, 1, 1, 2, 2]
    asubj = [y, z, x, y, z, y]
    aval  = [1.0, 1.0, 1.0, -1.0, 1.0, -1.0]
    bkc   = [boundkey.up, boundkey.ra, boundkey.ra]
    blc   = [-inf, log(alpha), log(gamma)]
    buc   = [log(Af), log(beta), log(delta)]

    with Task() as task:
        task.set_Stream(streamtype.log, streamprinter)

        # Add variables and constraints
        task.appendvars(numvar)
        task.appendcons(numcon)

        # Objective is the sum of three first variables
        task.putobjsense(objsense.maximize)
        task.putcslice(0, numvar, cval)
        task.putvarboundsliceconst(0, numvar, boundkey.fr, -inf, inf)

        # Add the three linear constraints
        task.putaijlist(asubi, asubj, aval)
        task.putconboundslice(0, numvar, bkc, blc, buc)

        # Affine expressions appearing in affine conic constraints
        # in this order:
        # u1, u2, x+y+log(2/Awall), x+z+log(2/Awall), 1.0, u1+u2-1.0
        numafe    = 6
        u1, u2 = 3, 4     # Indices of slack variables
        afeidx = [0, 1, 2, 2, 3, 3, 5, 5]
        varidx = [u1, u2, x, y, x, z, u1, u2]
        fval   = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
        gfull  = [0, 0, log(2/Aw), log(2/Aw), 1.0, -1.0]

        # New variables u1, u2
        task.appendvars(2)
        task.putvarboundsliceconst(u1, u2+1, boundkey.fr, -inf, inf)

        # Append affine expressions
        task.appendafes(numafe)
        task.putafefentrylist(afeidx, varidx, fval)
        task.putafegslice(0, numafe, gfull)

        # Two affine conic constraints
        expdom = task.appendprimalexpconedomain()

        # (u1, 1, x+y+log(2/Awall)) \in EXP
        task.appendacc(expdom, [0, 4, 2], None)

        # (u2, 1, x+z+log(2/Awall)) \in EXP
        task.appendacc(expdom, [1, 4, 3], None)

        # The constraint u1+u2-1 \in \ZERO is added also as an ACC
        task.appendacc(task.appendrzerodomain(1), [5], None)

        # Solve and map to original h, w, d
        task.optimize()
        task.writedata("gp1.ptf");
        xyz = task.getxxslice(soltype.itr, 0, numvar)
        return exp(xyz)

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

helloworld.py

Listing 17.11 helloworld.py Click here to download.
##
#  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File:      helloworld.py
#
#  The most basic example of how to get started with MOSEK.

from mosek import *

with Env() as env:                            # Create Environment
  with env.Task(0, 1) as task:                # Create Task
    task.appendvars(1)                          # 1 variable x
    task.putcj(0, 1.0)                          # c_0 = 1.0
    task.putvarbound(0, boundkey.ra, 2.0, 3.0)  # 2.0 <= x <= 3.0
    task.putobjsense(objsense.minimize)         # minimize

    task.optimize()                           # Optimize

    x = task.getxx(soltype.itr)               # Get solution
    print("Solution x = {}".format(x[0]))     # Print solution

lo1.py

Listing 17.12 lo1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      lo1.py
#
#  Purpose :   Demonstrates how to solve small linear
#              optimization problem using the MOSEK Python API.
##
import sys
import mosek

# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Create a task object
    with mosek.Task() as task:
        # Attach a log stream printer to the task
        task.set_Stream(mosek.streamtype.log, streamprinter)

        # Bound keys for constraints
        bkc = [mosek.boundkey.fx,
               mosek.boundkey.lo,
               mosek.boundkey.up]

        # Bound values for constraints
        blc = [30.0, 15.0, -inf]
        buc = [30.0, +inf, 25.0]

        # Bound keys for variables
        bkx = [mosek.boundkey.lo,
               mosek.boundkey.ra,
               mosek.boundkey.lo,
               mosek.boundkey.lo]

        # Bound values for variables
        blx = [0.0, 0.0, 0.0, 0.0]
        bux = [+inf, 10.0, +inf, +inf]

        # Objective coefficients
        c = [3.0, 1.0, 5.0, 1.0]

        # Below is the sparse representation of the A
        # matrix stored by column.
        asub = [[0, 1],
                [0, 1, 2],
                [0, 1],
                [1, 2]]
        aval = [[3.0, 2.0],
                [1.0, 1.0, 2.0],
                [2.0, 3.0],
                [1.0, 3.0]]

        numvar = len(bkx)
        numcon = len(bkc)

        # Append 'numcon' empty constraints.
        # The constraints will initially have no bounds.
        task.appendcons(numcon)

        # Append 'numvar' variables.
        # The variables will initially be fixed at zero (x=0).
        task.appendvars(numvar)

        for j in range(numvar):
            # Set the linear term c_j in the objective.
            task.putcj(j, c[j])

            # Set the bounds on variable j
            # blx[j] <= x_j <= bux[j]
            task.putvarbound(j, bkx[j], blx[j], bux[j])

            # Input column j of A
            task.putacol(j,                  # Variable (column) index.
                         asub[j],            # Row index of non-zeros in column j.
                         aval[j])            # Non-zero Values of column j.

        # Set the bounds on constraints.
         # blc[i] <= constraint_i <= buc[i]
        for i in range(numcon):
            task.putconbound(i, bkc[i], blc[i], buc[i])

        # Input the objective sense (minimize/maximize)
        task.putobjsense(mosek.objsense.maximize)

        # Solve the problem
        task.optimize()
        # Print a summary containing information
        # about the solution for debugging purposes
        task.solutionsummary(mosek.streamtype.msg)

        # Get status information about the solution
        solsta = task.getsolsta(mosek.soltype.bas)

        if (solsta == mosek.solsta.optimal):
            xx = task.getxx(mosek.soltype.bas)
            
            print("Optimal solution: ")
            for i in range(numvar):
                print("x[" + str(i) + "]=" + str(xx[i]))
        elif (solsta == mosek.solsta.dual_infeas_cer or
              solsta == mosek.solsta.prim_infeas_cer):
            print("Primal or dual infeasibility certificate found.\n")
        elif solsta == mosek.solsta.unknown:
            print("Unknown solution status")
        else:
            print("Other solution status")


# call the main function
try:
    main()
except mosek.Error as e:
    print("ERROR: %s" % str(e.errno))
    if e.msg is not None:
        print("\t%s" % e.msg)
        sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

lo2.py

Listing 17.13 lo2.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      lo2.py
#
#  Purpose :   Demonstrates how to solve small linear
#              optimization problem using the MOSEK Python API.
##
import sys
import mosek

# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Make a MOSEK environment
    with mosek.Env() as env:
        # Attach a printer to the environment
        env.set_Stream(mosek.streamtype.log, streamprinter)

        # Create a task
        with env.Task(0, 0) as task:
            # Attach a printer to the task
            task.set_Stream(mosek.streamtype.log, streamprinter)

            # Bound keys for constraints
            bkc = [mosek.boundkey.fx,
                   mosek.boundkey.lo,
                   mosek.boundkey.up]
            # Bound values for constraints
            blc = [30.0, 15.0, -inf]
            buc = [30.0, +inf, 25.0]
            # Bound keys for variables
            bkx = [mosek.boundkey.lo,
                   mosek.boundkey.ra,
                   mosek.boundkey.lo,
                   mosek.boundkey.lo]
            # Bound values for variables
            blx = [0.0, 0.0, 0.0, 0.0]
            bux = [+inf, 10.0, +inf, +inf]
            # Objective coefficients

            c = [3.0, 1.0, 5.0, 1.0]

            # We input the A matrix column-wise
            # asub contains row indexes
            asub = [[0, 1, 2],
                    [0, 1, 2, 3],
                    [0, 3]]
            # acof contains coefficients
            aval = [[3.0, 1.0, 2.0],
                    [2.0, 1.0, 3.0, 1.0],
                    [2.0, 3.0]]
            numvar = len(bkx)
            numcon = len(bkc)
            # Append 'numcon' empty constraints.
            # The constraints will initially have no bounds.
            task.appendcons(numcon)

            #Append 'numvar' variables.
            # The variables will initially be fixed at zero (x=0).
            task.appendvars(numvar)

            for j in range(numvar):
                # Set the linear term c_j in the objective.
                task.putcj(j, c[j])
                # Set the bounds on variable j
                # blx[j] <= x_j <= bux[j]
                task.putvarbound(j, bkx[j], blx[j], bux[j])

            for i in range(numcon):
                task.putconbound(i, bkc[i], blc[i], buc[i])
                # Input row i of A
                task.putarow(i,                     # Row index.
                             # Column indexes of non-zeros in row i.
                             asub[i],
                             aval[i])              # Non-zero Values of row i.

            # Input the objective sense (minimize/maximize)
            task.putobjsense(mosek.objsense.maximize)

            # Optimize the task
            task.optimize()

            # Print a summary containing information
            #about the solution for debugging purposes
            task.solutionsummary(mosek.streamtype.msg)

            prosta = task.getprosta(mosek.soltype.bas)
            solsta = task.getsolsta(mosek.soltype.bas)

            # Output a solution
            xx = task.getxx(mosek.soltype.bas)

            if solsta == mosek.solsta.optimal:
                print("Optimal solution: %s" % xx)
            elif solsta == mosek.solsta.dual_infeas_cer:
                print("Primal or dual infeasibility.\n")
            elif solsta == mosek.solsta.prim_infeas_cer:
                print("Primal or dual infeasibility.\n")
            elif mosek.solsta.unknown:
                print("Unknown solution status")
            else:
                print("Other solution status")

                

# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    if e.msg is not None:
        print("\t%s" % e.msg)
        sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

logistic.py

Listing 17.14 logistic.py Click here to download.
# Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File:      logistic.py
#
# Purpose: Implements logistic regression with regulatization.
#
#          Demonstrates using the exponential cone and log-sum-exp in Optimizer API.
#
#          Plots an example for 2D datasets.
from mosek import *
import numpy as np 
import sys, itertools

inf = 0.0

# Adds ACCs for t_i >= log ( 1 + exp((1-2*y[i]) * theta' * X[i]) )
# Adds auxiliary variables, AFE rows and constraints
def softplus(task, d, n, theta, t, X, y):
    nvar = task.getnumvar()
    ncon = task.getnumcon()
    nafe = task.getnumafe()
    task.appendvars(2*n)    # z1, z2
    task.appendcons(n)      # z1 + z2 = 1
    task.appendafes(4*n)    #theta * X[i] - t[i], -t[i], z1[i], z2[i]
    z1, z2 = nvar, nvar+n
    zcon = ncon
    thetaafe, tafe, z1afe, z2afe = nafe, nafe+n, nafe+2*n, nafe+3*n
    for i in range(n):
        task.putvarname(z1+i,f"z1[{i}]")
        task.putvarname(z2+i,f"z2[{i}]")

    # z1 + z2 = 1
    task.putaijlist(range(zcon, zcon+n), range(z1, z1+n), [1]*n)
    task.putaijlist(range(zcon, zcon+n), range(z2, z2+n), [1]*n)
    task.putconboundsliceconst(zcon, zcon+n, boundkey.fx, 1, 1)
    task.putvarboundsliceconst(nvar, nvar+2*n, boundkey.fr, -inf, inf)

    # Affine conic expressions
    afeidx, varidx, fval = [], [], []

    ## Thetas
    for i in range(n):
      for j in range(d):
        afeidx.append(thetaafe + i)
        varidx.append(theta + j)
        fval.append(-X[i][j] if y[i]==1 else X[i][j])

    # -t[i]
    afeidx.extend([thetaafe + i for i in range(n)] + [tafe + i for i in range(n)])
    varidx.extend([t + i for i in range(n)] + [t + i for i in range(n)])
    fval.extend([-1.0]*(2*n))

    # z1, z2
    afeidx.extend([z1afe + i for i in range(n)] + [z2afe + i for i in range(n)])
    varidx.extend([z1 + i for i in range(n)] + [z2 + i for i in range(n)])
    fval.extend([1.0]*(2*n))

    # Add the expressions
    task.putafefentrylist(afeidx, varidx, fval)

    # Add a single row with the constant expression "1.0"
    oneafe = task.getnumafe()
    task.appendafes(1)
    task.putafeg(oneafe, 1.0)

    # Add an exponential cone domain
    expdomain = task.appendprimalexpconedomain()

    # Conic constraints
    acci = task.getnumacc()
    for i in range(n):
      task.appendacc(expdomain, [z1afe+i, oneafe, thetaafe+i], None)
      task.appendacc(expdomain, [z2afe+i, oneafe, tafe+i], None)
      task.putaccname(acci,  f"z1:theta[{i}]")
      task.putaccname(acci+1,f"z2:t[{i}]")
      acci += 2

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

    with env.Task() as task:
        # Variables [r; theta; t; u]
        nvar = 1+d+2*n
        task.appendvars(nvar)
        task.putvarboundsliceconst(0, nvar, boundkey.fr, -inf, inf)
        r, theta, t = 0, 1, 1+d
        task.putvarname(r,"r");
        for j in range(d): task.putvarname(theta+j,f"theta[{j}]");
        for j in range(n): task.putvarname(t+j,f"t[{j}]");

        # Objective lambda*r + sum(t)
        task.putobjsense(objsense.minimize)
        task.putcj(r, lamb)
        task.putclist(range(t, t+n), [1.0]*n)

        # Softplus function constraints
        softplus(task, d, n, theta, t, X, y);

        # Regularization
        # Append a sequence of linear expressions (r, theta) to F
        numafe = task.getnumafe()
        task.appendafes(1+d)
        task.putafefentry(numafe, r, 1.0)
        for i in range(d):
          task.putafefentry(numafe + i + 1, theta + i, 1.0)

        # Add the constraint
        task.appendaccseq(task.appendquadraticconedomain(1+d), numafe, None)

        # Solution
        task.writedata('logistic.ptf')
        task.optimize()
        xx = task.getxxslice(soltype.itr, theta, theta+d)

        return xx

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

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

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

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

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

    plt.show()

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

env = Env()

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

# Example 2: discover a circle
x, y = [], []
for x1,x2 in itertools.product(np.linspace(-1, 1, 30),np.linspace(-1, 1, 30)):
    x.append([x1,x2])
    y.append(0 if x1**2+x2**2<0.69 else 1)

x, y = np.array(x), np.array(y)
X = mapFeatures(x, d=2)
# print(list(x[:,0]))
# print(list(x[:,1]))
# for i in range(len(x)):
#     print(f"X[{i}] = {list(X[i,:])}")

theta = logisticRegression(env, X, y, 0.1)
print(theta)

mico1.py

Listing 17.15 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 mosek, sys

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()

with mosek.Task() as task:
    task.set_Stream(mosek.streamtype.log, streamprinter)

    task.appendvars(3)    # x, y, t
    x, y, t = 0, 1, 2
    task.putvarboundsliceconst(0, 3, mosek.boundkey.fr, -0.0, 0.0)
    
    # Integrality constraints
    task.putvartypelist([x,y], [mosek.variabletype.type_int]*2)

    # Set up the affine expression
    # x, x-3.8, y, t, 1.0
    task.appendafes(5)
    task.putafefentrylist([0, 1, 2, 3],
                          [x,x,y,t],
                          [1,1,1,1])
    task.putafegslice(0, 5, [0, -3.8, 0, 0, 1.0])

    # Add constraint (x-3.8, 1, y) \in \EXP
    task.appendacc(task.appendprimalexpconedomain(), [1, 4, 2], None)

    # Add constraint (t, x, y) \in \QUAD
    task.appendacc(task.appendquadraticconedomain(3), [3, 0, 2], None)

    # Objective
    task.putobjsense(mosek.objsense.minimize)
    task.putcj(t, 1)

    # Optimize the task
    task.optimize()
    task.solutionsummary(mosek.streamtype.msg)

    xx = task.getxxslice(mosek.soltype.itg, 0, 2)
    print(xx)

milo1.py

Listing 17.16 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 using the MOSEK Python API.
##
import sys
import mosek

# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Make a MOSEK environment
    with mosek.Env() as env:
        # Attach a printer to the environment
        env.set_Stream(mosek.streamtype.log, streamprinter)

        # Create a task
        with env.Task(0, 0) as task:
            # Attach a printer to the task
            task.set_Stream(mosek.streamtype.log, streamprinter)

            bkc = [mosek.boundkey.up, mosek.boundkey.lo]
            blc = [-inf, -4.0]
            buc = [250.0, inf]

            bkx = [mosek.boundkey.lo, mosek.boundkey.lo]
            blx = [0.0, 0.0]
            bux = [inf, inf]

            c = [1.0, 0.64]

            asub = [[0, 1], [0, 1]]
            aval = [[50.0, 3.0], [31.0, -2.0]]

            numvar = len(bkx)
            numcon = len(bkc)

            # Append 'numcon' empty constraints.
            # The constraints will initially have no bounds.
            task.appendcons(numcon)

            #Append 'numvar' variables.
            # The variables will initially be fixed at zero (x=0).
            task.appendvars(numvar)

            for j in range(numvar):
                # Set the linear term c_j in the objective.
                task.putcj(j, c[j])
                # Set the bounds on variable j
                # blx[j] <= x_j <= bux[j]
                task.putvarbound(j, bkx[j], blx[j], bux[j])
                # Input column j of A
                task.putacol(j,                  # Variable (column) index.
                             # Row index of non-zeros in column j.
                             asub[j],
                             aval[j])            # Non-zero Values of column j.

            task.putconboundlist(range(numcon), bkc, blc, buc)

            # Input the objective sense (minimize/maximize)
            task.putobjsense(mosek.objsense.maximize)

            # Define variables to be integers
            task.putvartypelist([0, 1],
                                [mosek.variabletype.type_int,
                                 mosek.variabletype.type_int])

            # Set max solution time
            task.putdouparam(mosek.dparam.mio_max_time, 60.0);

            # Optimize the task
            task.optimize()
            task.writedata("milo1.ptf")

            # Print a summary containing information
            # about the solution for debugging purposes
            task.solutionsummary(mosek.streamtype.msg)

            prosta = task.getprosta(mosek.soltype.itg)
            solsta = task.getsolsta(mosek.soltype.itg)

            # Output a solution         
            xx = task.getxx(mosek.soltype.itg)

            if solsta in [mosek.solsta.integer_optimal]:
                print("Optimal solution: %s" % xx)
            elif solsta == mosek.solsta.prim_feas:
                print("Feasible solution: %s" % xx)
            elif mosek.solsta.unknown:
                if prosta == mosek.prosta.prim_infeas_or_unbounded:
                    print("Problem status Infeasible or unbounded.\n")
                elif prosta == mosek.prosta.prim_infeas:
                    print("Problem status Infeasible.\n")
                elif prosta == mosek.prosta.unkown:
                    print("Problem status unkown.\n")
                else:
                    print("Other problem status.\n")
            else:
                print("Other solution status")


# call the main function
try:
    main()
except mosek.MosekException as msg:
    #print "ERROR: %s" % str(code)
    if msg is not None:
        print("\t%s" % msg)
        sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

mioinitsol.py

Listing 17.17 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 using the MOSEK Python API.
##
import sys
import mosek

# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Make a MOSEK environment
    with mosek.Env() as env:
        # Attach a printer to the environment
        env.set_Stream(mosek.streamtype.log, streamprinter)

        # Create a task
        with env.Task(0, 0) as task:
            # Attach a printer to the task
            task.set_Stream(mosek.streamtype.log, streamprinter)

            bkc = [mosek.boundkey.up]
            blc = [-inf, ]
            buc = [2.5]

            bkx = [mosek.boundkey.lo,
                   mosek.boundkey.lo,
                   mosek.boundkey.lo,
                   mosek.boundkey.lo]

            blx = [0.0, 0.0, 0.0, 0.0]
            bux = [inf, inf, inf, inf]

            c = [7.0, 10.0, 1.0, 5.0]

            asub = [0, 0, 0, 0]
            acof = [1.0, 1.0, 1.0, 1.0]

            ptrb = [0, 1, 2, 3]
            ptre = [1, 2, 3, 4]

            numvar = len(bkx)
            numcon = len(bkc)

            # Input linear data
            task.inputdata(numcon, numvar,
                           c, 0.0,
                           ptrb, ptre, asub, acof,
                           bkc, blc, buc,
                           bkx, blx, bux)

            # Input objective sense
            task.putobjsense(mosek.objsense.maximize)

            # Define variables to be integers
            task.putvartypelist([0, 1, 2],
                                [mosek.variabletype.type_int]*3)

            # Assign values to integer variables. 
            # (We only set a slice of xx)
            task.putxxslice(mosek.soltype.itg, 0, 3, [1.0, 1.0, 0.0])

            # Request constructing the solution from integer variable values
            task.putintparam(mosek.iparam.mio_construct_sol, mosek.onoffkey.on)

            # Optimize
            task.optimize()
            task.writedata("mioinitsol.ptf")

            task.solutionsummary(mosek.streamtype.log)

            if task.solutiondef(mosek.soltype.itg):
                # Output a solution
                xx = task.getxx(mosek.soltype.itg)

                print("Integer optimal solution")
                for j in range(0, numvar):
                    print("\tx[%d] = %e" % (j, xx[j]))

                # Was the initial guess used?
                constr = task.getintinf(mosek.iinfitem.mio_construct_solution)
                constrVal = task.getdouinf(mosek.dinfitem.mio_construct_solution_obj)
                print("Construct solution utilization: {0}\nConstruct solution objective: {1:.3f}\n".format(constr, constrVal))
            else:
                print("No integer solution is available.")


# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    if e.msg is not None:
        print("\t%s" % e.msg)
    sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

opt_server_async.py

Listing 17.18 opt_server_async.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      opt_server_async.py
#
#  Purpose :   Demonstrates how to use MOSEK OptServer
#              to solve optimization problem asynchronously
##
import mosek
import sys
import time

def streamprinter(msg):
    sys.stdout.write(msg)
    sys.stdout.flush()

if len(sys.argv) < 4:
    print("Missing argument, syntax is:")
    print("  opt-server-async inputfile host:port numpolls [cert]")
else:
    filename = sys.argv[1]
    addr = sys.argv[2]
    numpolls = int(sys.argv[3])
    token = None
    cert = None if len(sys.argv) < 5 else sys.argv[4]


    with mosek.Env() as env:

        with env.Task(0, 0) as task:

            print("reading task from file")
            task.readdata(filename)

            if cert is not None:
                task.putstrparam(mosek.sparam.remote_tls_cert_path,cert)

            print("Solve the problem remotely (async)")
            token = task.asyncoptimize(addr,"")

        print("Task token: %s" % token)

        with env.Task(0, 0) as task:

            task.readdata(filename)

            if cert is not None:
                task.putstrparam(mosek.sparam.remote_tls_cert_path,cert)
            task.set_Stream(mosek.streamtype.log, streamprinter)

            i = 0

            while i < numpolls:

                time.sleep(0.1)

                print("poll %d..." % i)
                respavailable, res, trm = task.asyncpoll(addr,
                                                         "",
                                                         token)

                print("done!")

                if respavailable:
                    print("solution available!")

                    respavailable, res, trm = task.asyncgetresult(addr,
                                                                  "",
                                                                  token)

                    task.solutionsummary(mosek.streamtype.log)
                    break

                i = i + 1

                if i == numpolls:
                    print("max number of polls reached, stopping host.")
                    task.asyncstop(addr,"", token)

opt_server_sync.py

Listing 17.19 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
##
import mosek
import sys

def streamprinter(msg):
    sys.stdout.write(msg)
    sys.stdout.flush()

if len(sys.argv) <= 2:
    print("Missing argument, syntax is:")
    print("  opt_server_sync inputfile addr [certpath]")
else:

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

    # Create the mosek environment.
    with mosek.Env() as env:

        # Create a task object linked with the environment env.
        # We create it with 0 variables and 0 constraints initially,
        # since we do not know the size of the problem.
        with env.Task(0, 0) as task:
            task.set_Stream(mosek.streamtype.log, streamprinter)

            # We assume that a problem file was given as the first command
            # line argument (received in `argv')
            task.readdata(inputfile)
            
            # Set OptServer URL
            task.putoptserverhost(serveraddr)

            # Path to certificate, if any
            if tlscert is not None:
                task.putstrparam(mosek.sparam.remote_tls_cert_path, tlscert)

            # Solve the problem remotely, no access token
            trm = task.optimize()

            # Print a summary of the solution
            task.solutionsummary(mosek.streamtype.log)

parallel.py

Listing 17.20 parallel.py Click here to download.
#
#  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File:      parallel.py
#
#  Purpose: Demonstrates parallel optimization using optimizebatch()
#
import mosek, sys
 
# Example of how to use env.optimizebatch()
# Optimizes tasks whose names were read from command line.
def main(argv):
  n = len(argv) - 1
  tasks = []

  threadpoolsize = 6                   # Size of thread pool available for all tasks

  with mosek.Env() as env:
    # Set up some example list of tasks to optimize
    for i in range(n):
      t = mosek.Task(env, 0, 0)
      t.readdata(argv[i+1])
      
      # We can set the number of threads for each task
      t.putintparam(mosek.iparam.num_threads, 2)
      tasks.append(t)

    # Optimize all the given tasks in parallel
    trm, res = env.optimizebatch(False,               # No race
                                 -1.0,                # No time limit
                                 threadpoolsize,
                                 tasks)               # List of tasks to optimize

    for i in range(n):
      print("Task  {0}  res {1}   trm {2}   obj_val  {3}  time {4}".format(
             i, 
             res[i], 
             trm[i],
             tasks[i].getdouinf(mosek.dinfitem.intpnt_primal_obj),
             tasks[i].getdouinf(mosek.dinfitem.optimizer_time)))

if __name__ == "__main__":
    main(sys.argv)

parameters.py

Listing 17.21 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 get/set
#              parameters with MOSEK Python API
##
import mosek

# Create the mosek environment.
with mosek.Env() as env:
    # Create the mosek environment.
    with env.Task() as task:
        print('Test MOSEK parameter get/set functions')

        # Set log level (integer parameter)
        task.putintparam(mosek.iparam.log, 1)
        # Select interior-point optimizer... (integer parameter)
        task.putintparam(mosek.iparam.optimizer, mosek.optimizertype.intpnt)
        # ... without basis identification (integer parameter)
        task.putintparam(mosek.iparam.intpnt_basis, mosek.basindtype.never)
        # Set relative gap tolerance (double parameter)
        task.putdouparam(mosek.dparam.intpnt_co_tol_rel_gap, 1.0e-7)

        # The same using explicit string names 
        task.putparam     ("MSK_DPAR_INTPNT_CO_TOL_REL_GAP", "1.0e-7")
        task.putnadouparam("MSK_DPAR_INTPNT_CO_TOL_REL_GAP",  1.0e-7 )

        # Incorrect value
        try:
            task.putdouparam(mosek.dparam.intpnt_co_tol_rel_gap, -1.0)
        except:
            print('Wrong parameter value') 


        param = task.getdouparam(mosek.dparam.intpnt_co_tol_rel_gap)
        print('Current value for parameter intpnt_co_tol_rel_gap = {}'.format(param))

        # Define and solve an optimization problem here
        # task.optimize() 
        # After optimization: 

        print('Get MOSEK information items')

        tm = task.getdouinf(mosek.dinfitem.optimizer_time)
        it = task.getintinf(mosek.iinfitem.intpnt_iter)

        print('Time: {0}\nIterations: {1}'.format(tm,it))

pinfeas.py

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

# Set up a simple linear problem from the manual for test purposes
def testProblem():
    inf = 0.0
    task = Task()
    task.appendvars(7)
    task.appendcons(7)
    task.putclist(range(0, 7), [1, 2, 5, 2, 1, 2, 1])
    task.putaijlist([0,0,1,1,2,2,2,3,3,4,5,5,6,6],
                    [0,1,2,3,4,5,6,0,4,1,2,5,3,6],
                    [1] * 14)
    task.putconboundslice(0, 7, [boundkey.up]*3+[boundkey.fx]*4,
                                [-inf, -inf, -inf, 1100, 200, 500, 500],
                                [200, 1000, 1000, 1100, 200, 500, 500])
    task.putvarboundsliceconst(0, 7, boundkey.lo, 0, +inf)
    return task

# Analyzes and prints infeasibility contributing elements
# sl - dual values for lower bounds
# su - dual values for upper bounds
# eps - tolerance for when a nunzero dual value is significant
def analyzeCertificate(sl, su, eps):
    n = len(sl)
    for i in range(n):
        if abs(sl[i]) > eps:
            print(f"#{i}: lower, dual = {sl[i]}") 
        if abs(su[i]) > eps: 
            print(f"#{i}: upper, dual = {su[i]}") 

def pinfeas():
    # In this example we set up a simple problem
    # One could use any task or a task read from a file
    task = testProblem()

    # Useful for debugging
    task.writedata("pinfeas.ptf")
    task.set_Stream(streamtype.log, sys.stdout.write) 
    
    # Perform the optimization.
    task.optimize()
    task.solutionsummary(streamtype.log)

    # Check problem status, we use the interior point solution
    if task.getprosta(soltype.itr) == prosta.prim_infeas:
        # Set the tolerance at which we consider a dual value as essential
        eps = 1e-7

        print("Variable bounds important for infeasibility: ")
        analyzeCertificate(task.getslx(soltype.itr), task.getsux(soltype.itr), eps)
        
        print("Constraint bounds important for infeasibility: ")
        analyzeCertificate(task.getslc(soltype.itr), task.getsuc(soltype.itr), eps)
    else:
        print("The problem is not primal infeasible, no certificate to show");

pinfeas()

portfolio_1_basic.py

Listing 17.23 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.
##
import mosek
import sys
import numpy as np

if __name__ == '__main__':

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

    inf = 0.0 # This value has no significance

    # Variable offsets
    numvar = n
    voff_x = 0

    # Constraints offsets
    numcon = 1
    coff_bud = 0

    with mosek.Env() as env:
        with env.Task(0, 0) as task:
            task.set_Stream(mosek.streamtype.log, sys.stdout.write)

            # Holding variable x of length n
            # No other auxiliary variables are needed in this formulation
            task.appendvars(numvar)

            # Optionally we can give the variables names
            for j in range(0, n):
                task.putvarname(voff_x + j, "x[%d]" % (1 + j))

            # No short-selling in this model, all of x >= 0
            task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)

            # One linear constraint: total budget
            task.appendcons(1)
            task.putconname(coff_bud, "budget")
            task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n)        # e^T x
            rtemp = w + sum(x0)
            task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp)    # equals w + sum(x0)

            # Input (gamma, GTx) in the AFE (affine expression) storage
            # We need k+1 rows
            task.appendafes(k + 1)
            # The first affine expression = gamma
            task.putafeg(0, gamma)
            # The remaining k expressions comprise GT*x, we add them row by row
            # In more realisic scenarios it would be better to extract nonzeros and input in sparse form
            for i in range(0, k):
                task.putafefrow(i + 1, range(voff_x, voff_x + n), GT[i])

            # Input the affine conic constraint (gamma, GT*x) \in QCone
            # Add the quadratic domain of dimension k+1
            qdom = task.appendquadraticconedomain(k + 1)
            # Add the constraint
            task.appendaccseq(qdom, 0, None)
            task.putaccname(0, "risk")

            # Objective: maximize expected return mu^T x
            task.putclist(range(voff_x, voff_x + n), mu)
            task.putobjsense(mosek.objsense.maximize)

            # Dump the problem to a human readable PTF file.
            task.writedata("dump.ptf")

            # Solve the problem
            task.optimize()

            # Display solution summary for quick inspection of results.
            # In this simplified example we skip checks for problem and solution status
            task.solutionsummary(mosek.streamtype.msg)

            # Check if the interior point solution is an optimal point
            solsta = task.getsolsta(mosek.soltype.itr)
            if (solsta != mosek.solsta.optimal):
                # See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
                raise Exception(f"Unexpected solution status: {solsta}")

            # Retrieve results
            xx     = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
            expret = task.getprimalobj(mosek.soltype.itr)

            print(f'Expected return: {expret:.10e} Std. deviation: {gamma:.4e}')
            np.set_printoptions(precision=4)
            print(f'Optimal portfolio: {np.array(xx)}')

portfolio_2_frontier.py

Listing 17.24 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.
##
import mosek
import sys
import numpy as np

if __name__ == '__main__':

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

    alphas = [0.0, 0.01, 0.1, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 10.0]

    inf = 0.0 # This value has no significance

    # Offset of variables into the API variable.
    numvar = n + 1
    voff_x = 0
    voff_s = n

    # Offset of constraints
    coff_bud = 0

    with mosek.Env() as env:
        with env.Task(0, 0) as task:
            task.set_Stream(mosek.streamtype.log, sys.stdout.write)

            # Variables: 
            task.appendvars(numvar)

            # Optionally we can give the variables names
            for j in range(0, n):
                task.putvarname(voff_x + j, "x[%d]" % (1 + j))
            task.putvarname(voff_s, "s")

            # No short-selling in this model, all of x >= 0
            task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)

            # s is free variable
            task.putvarbound(voff_s, mosek.boundkey.fr, -inf, inf)

            # One linear constraint: total budget
            task.appendcons(1)
            task.putconname(coff_bud, "budget")
            task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n)        # e^T x
            rtemp = w + sum(x0)
            task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp)    # equals w + sum(x0)

            # Input (gamma, GTx) in the AFE (affine expression) storage
            # We build the following F and g for variables [x, s]:
            #     [0, 1]      [0  ]
            # F = [0, 0], g = [0.5]
            #     [GT,0]      [0  ]
            # We need k+2 rows
            task.appendafes(k + 2)
            # The first affine expression is variable s (last variable, index n)
            task.putafefrow(0, [n], [1.0])            
            # The second affine expression is constant 0.5
            task.putafeg(1, 0.5)
            # The remaining k affine expressions comprise GT*x, we add them row by row
            # In more realisic scenarios it would be better to extract nonzeros and input in sparse form
            for i in range(0, k):
                task.putafefrow(i + 2, range(voff_x, voff_x + n), GT[i])

            # Input the affine conic constraint (s, 0.5, GT*x) \in RQCone
            # Add the quadratic domain of dimension k+1
            rqdom = task.appendrquadraticconedomain(k + 2)
            # Add the constraint
            task.appendaccseq(rqdom, 0, None)
            task.putaccname(0, "risk")

            # Set objective coefficients (x part): mu'x - alpha * s
            task.putclist(range(voff_x, voff_x + n), mu)

            task.putobjsense(mosek.objsense.maximize)

            # Turn all log output off.
            task.putintparam(mosek.iparam.log, 0)

            for alpha in alphas:
                # Dump the problem to a human readable PTF file.
                task.writedata("dump.ptf")

                task.putcj(voff_s, -alpha)

                task.optimize()

                # Display the solution summary for quick inspection of results.
                # task.solutionsummary(mosek.streamtype.msg)

                # Check if the interior point solution is an optimal point
                solsta = task.getsolsta(mosek.soltype.itr)

                if solsta in [mosek.solsta.optimal]:
                    expret = 0.0
                    x = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
                    for j in range(0, n):
                        expret += mu[j] * x[j]

                    stddev = np.sqrt(task.getxxslice(mosek.soltype.itr, voff_s, voff_s + 1))

                    print("alpha = {0:.2e} exp. ret. = {1:.3e}, std. dev. {2:.3e}".format(alpha, expret, stddev[0]))
                else:
                    # See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
                    print("An error occurred when solving for alpha=%e\n" % alpha)
                    raise Exception(f"Unexpected solution status: {solsta}")

portfolio_3_impact.py

Listing 17.25 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.
##
import mosek
import sys

if __name__ == '__main__':

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

    # This value has no significance.
    inf = 0.0

    # Offset of variables.
    numvar = 3 * n
    voff_x, voff_c, voff_z = 0, n, 2 * n

    # Offset of constraints.
    numcon = 2 * n + 1
    coff_bud, coff_abs1, coff_abs2 = 0, 1, 1 + n

    with mosek.Env() as env:
        with env.Task(0, 0) as task:
            task.set_Stream(mosek.streamtype.log, sys.stdout.write)

            # Variables (vector of x, c, z)
            task.appendvars(numvar)
            for j in range(0, n):
                task.putvarname(voff_x + j, "x[%d]" % (j + 1))
                task.putvarname(voff_c + j, "c[%d]" % (j + 1))
                task.putvarname(voff_z + j, "z[%d]" % (j + 1))
            
            # Apply variable bounds (x >= 0, c and z free)
            task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
            task.putvarboundsliceconst(voff_c, voff_c + n, mosek.boundkey.fr, -inf, inf)
            task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)

            # Linear constraints
            # - Budget
            task.appendcons(1)
            task.putconname(coff_bud, "budget")
            task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n)    # e^T x
            task.putaijlist([coff_bud] * n, range(voff_c, voff_c + n), m)            # m^T c
            rtemp = w + sum(x0)
            task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp)    # equals w + sum(x0)

            # - Absolute value
            task.appendcons(2 * n)
            for i in range(0, n):
                task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
                task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
            task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), n * [-1.0])
            task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), n * [1.0])
            task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)         
            task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), n * [1.0])
            task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), n * [1.0])
            task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)            

            # ACCs
            aoff_q = 0
            aoff_pow = k + 1
            # - (gamma, GTx) in Q(k+1)
            # The part of F and g for variable x:
            #     [0,  0, 0]      [gamma]
            # F = [GT, 0, 0], g = [0    ]    
            task.appendafes(k + 1)
            task.putafeg(aoff_q, gamma)
            for i in range(0, k):
                task.putafefrow(aoff_q + 1 + i, range(voff_x, voff_x + n), GT[i])
            qdom = task.appendquadraticconedomain(k + 1)
            task.appendaccseq(qdom, aoff_q, None)
            task.putaccname(aoff_q, "risk")

            # - (c_j, 1, z_j) in P3(2/3, 1/3)
            # The part of F and g for variables [c, z]:
            #     [0, I, 0]      [0]
            # F = [0, 0, I], g = [0]
            #     [0, 0, 0]      [1]
            task.appendafes(2 * n + 1)
            task.putafefentrylist(range(aoff_pow, aoff_pow + n), range(voff_c, voff_c + n), [1.0] * n)
            task.putafefentrylist(range(aoff_pow + n, aoff_pow + 2 * n), range(voff_z, voff_z + n), [1.0] * n)
            task.putafeg(aoff_pow + 2 * n, 1.0)
            # We use one row from F and g for both c_j and z_j, and the last row of F and g for the constant 1.
            # NOTE: Here we reuse the last AFE and the power cone n times, but we store them only once.
            powdom = task.appendprimalpowerconedomain(3, [2, 1])
            afe_list = [(aoff_pow + i, aoff_pow + 2 * n, aoff_pow + n + i) for i in range(0, n)]
            flat_afe_list = [idx for sublist in afe_list for idx in sublist]
            task.appendaccs([powdom] * n, flat_afe_list, None)
            for i in range(0, n): 
                task.putaccname(i + 1, "market_impact[%d]" % i)
            
            # Objective
            task.putclist(range(voff_x, voff_x + n), mu)      
            task.putobjsense(mosek.objsense.maximize)

            # Turn all log output off.
            # task.putintparam(mosek.iparam.log,0)

            # Dump the problem to a human readable PTF file.
            task.writedata("dump.ptf")

            task.optimize()

            # Display the solution summary for quick inspection of results.
            task.solutionsummary(mosek.streamtype.msg)

            # Check if the interior point solution is an optimal point
            solsta = task.getsolsta(mosek.soltype.itr)
            if (solsta != mosek.solsta.optimal):
                # See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
                raise Exception(f"Unexpected solution status: {solsta}")

            expret = 0.0
            x = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
            for j in range(0, n):
                expret += mu[j] * x[j]

            print("\nExpected return %e for gamma %e\n" % (expret, gamma))

portfolio_4_transcost.py

Listing 17.26 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.
##
import mosek
import sys

if __name__ == '__main__':

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

    f = n*[0.01]
    g = n*[0.001]
    gamma = 0.36

    # This value has no significance.
    inf = 0.0

    # Offset of variables.
    numvar = 3 * n
    voff_x, voff_z, voff_y = 0, n, 2 * n

    # Offset of constraints.
    numcon = 3 * n + 1
    coff_bud, coff_abs1, coff_abs2, coff_swi = 0, 1, 1 + n, 1 + 2 * n 

    with mosek.Env() as env:
        with env.Task(0, 0) as task:
            task.set_Stream(mosek.streamtype.log, sys.stdout.write)

            # Variables (vector of x, z, y)
            task.appendvars(numvar)
            for j in range(0, n):
                task.putvarname(voff_x + j, "x[%d]" % (j + 1))
                task.putvarname(voff_z + j, "z[%d]" % (j + 1))
                task.putvarname(voff_y + j, "y[%d]" % (j + 1))
            
            # Apply variable bounds (x >= 0, z free, y binary)
            task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
            task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)
            task.putvarboundsliceconst(voff_y, voff_y + n, mosek.boundkey.ra, 0.0, 1.0)
            task.putvartypelist(range(voff_y, voff_y + n), [mosek.variabletype.type_int] * n)

            # Linear constraints
            # - Budget
            task.appendcons(1)
            task.putconname(coff_bud, "budget")
            task.putaijlist([0] * n, range(voff_x, voff_x + n), [1.0] * n)    # e^T x
            task.putaijlist([0] * n, range(voff_z, voff_z + n), g)            # g^T z
            task.putaijlist([0] * n, range(voff_y, voff_y + n), f)            # f^T y
            U = w + sum(x0)
            task.putconbound(coff_bud, mosek.boundkey.fx, U, U)    # equals w + sum(x0)

            # - Absolute value
            task.appendcons(2 * n)
            for i in range(0, n):
                task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
                task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
            task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), [-1.0] * n)
            task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), [1.0] * n)
            task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)         
            task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), [1.0] * n)
            task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), [1.0] * n)
            task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)      

            # - Switch 
            task.appendcons(n)
            for i in range(0, n):
                task.putconname(coff_swi + i, "switch[%d]" % (1 + i))
            task.putaijlist(range(coff_swi, coff_swi + n), range(voff_z, voff_z + n), [1.0] * n)         
            task.putaijlist(range(coff_swi, coff_swi + n), range(voff_y, voff_y + n), [-U] * n)
            task.putconboundlist(range(coff_swi, coff_swi + n), [mosek.boundkey.up] * n, [-inf] * n, [0.0] * n)      

            # ACCs
            aoff_q = 0
            # - (gamma, GTx) in Q(k+1)
            # The part of F and g for variable x:
            #     [0,  0, 0]      [gamma]
            # F = [GT, 0, 0], g = [0    ]    
            task.appendafes(k + 1)
            task.putafeg(aoff_q, gamma)
            for i in range(0, k):
                task.putafefrow(aoff_q + i + 1, range(voff_x, voff_x + n), GT[i])
            qdom = task.appendquadraticconedomain(k + 1)
            task.appendaccseq(qdom, aoff_q, None)
            task.putaccname(0, "risk")

            # Objective
            task.putclist(range(voff_x, voff_x + n), mu)      
            task.putobjsense(mosek.objsense.maximize)

            # Turn all log output off.
            # task.putintparam(mosek.iparam.log,0)

            # Dump the problem to a human readable PTF file.
            task.writedata("dump.ptf")

            task.optimize()

            # Display the solution summary for quick inspection of results.
            task.solutionsummary(mosek.streamtype.msg)

            # Check if the integer solution is an optimal point
            solsta = task.getsolsta(mosek.soltype.itg)
            if (solsta != mosek.solsta.integer_optimal):
                # See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
                raise Exception(f"Unexpected solution status: {solsta}")

            expret = 0.0
            x = task.getxxslice(mosek.soltype.itg, voff_x, voff_x + n)
            for j in range(0, n):
                expret += mu[j] * x[j]

            tcost = 0.0
            z = task.getxxslice(mosek.soltype.itg, voff_z, voff_z + n)
            y = task.getxxslice(mosek.soltype.itg, voff_y, voff_y + n)
            for j in range(0, n):
                tcost += g[j] * z[j] + f[j] * y[j]

            print("\nExpected return %e for gamma %e. Transaction cost: %e\n" % (expret, gamma, tcost))

portfolio_5_card.py

Listing 17.27 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.
##
import mosek
import sys
import numpy as np

# This value has no significance.
inf = 0.0 

def markowitz_with_card(n, k, x0, w, gamma, mu, GT, K):
    with mosek.Env() as env:
        with env.Task(0, 0) as task:
            task.set_Stream(mosek.streamtype.log, sys.stdout.write)

            # Offset of variables.
            numvar = 3 * n
            voff_x, voff_z, voff_y = 0, n, 2 * n

            # Offset of constraints.
            numcon = 3 * n + 2
            coff_bud, coff_abs1, coff_abs2, coff_swi, coff_card = 0, 1, 1 + n, 1 + 2 * n, 1 + 3 * n 

            # Variables (vector of x, z, y)
            task.appendvars(numvar)
            for j in range(0, n):
                task.putvarname(voff_x + j, "x[%d]" % (j + 1))
                task.putvarname(voff_z + j, "z[%d]" % (j + 1))
                task.putvarname(voff_y + j, "y[%d]" % (j + 1))
            
            # Apply variable bounds (x >= 0, z free, y binary)
            task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
            task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)
            task.putvarboundsliceconst(voff_y, voff_y + n, mosek.boundkey.ra, 0.0, 1.0)
            task.putvartypelist(range(voff_y, voff_y + n), [mosek.variabletype.type_int] * n)

            # Linear constraints
            # - Budget
            task.appendcons(1)
            task.putconname(coff_bud, "budget")
            task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n)    # e^T x
            U = w + sum(x0)
            task.putconbound(coff_bud, mosek.boundkey.fx, U, U)    # = w + sum(x0)

            # - Absolute value
            task.appendcons(2 * n)
            for i in range(0, n):
                task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
                task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
            task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), [-1.0] * n)
            task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), [1.0] * n)
            task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)         
            task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), [1.0] * n)
            task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), [1.0] * n)
            task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)      

            # - Switch 
            task.appendcons(n)
            for i in range(0, n):
                task.putconname(coff_swi + i, "switch[%d]" % (1 + i))
            task.putaijlist(range(coff_swi, coff_swi + n), range(voff_z, voff_z + n), [1.0] * n)         
            task.putaijlist(range(coff_swi, coff_swi + n), range(voff_y, voff_y + n), [-U] * n)
            task.putconboundlist(range(coff_swi, coff_swi + n), [mosek.boundkey.up] * n, [-inf] * n, [0.0] * n)      

            # - Cardinality
            task.appendcons(1)
            task.putconname(coff_card, "cardinality")
            task.putaijlist([coff_card] * n, range(voff_y, voff_y + n), [1.0] * n)    # e^T y
            task.putconbound(coff_card, mosek.boundkey.up, -inf, K)           # <= K

            # ACCs
            aoff_q = 0
            # - (gamma, GTx) in Q(k+1)
            # The part of F and g for variable x:
            #     [0,  0, 0]      [gamma]
            # F = [GT, 0, 0], g = [0    ]    
            task.appendafes(k + 1)
            task.putafeg(aoff_q, gamma)
            for i in range(0, k):
                task.putafefrow(aoff_q + i + 1, range(voff_x, voff_x + n), GT[i])
            qdom = task.appendquadraticconedomain(k + 1)
            task.appendaccseq(qdom, aoff_q, None)
            task.putaccname(0, "risk")

            # Objective
            task.putclist(range(voff_x, voff_x + n), mu)      
            task.putobjsense(mosek.objsense.maximize)

            # Turn all log output off.
            task.putintparam(mosek.iparam.log,0)

            # Dump the problem to a human readable PTF file.
            task.writedata("dump.ptf")

            task.optimize()

            # Check if the integer solution is an optimal point
            solsta = task.getsolsta(mosek.soltype.itg)
            if (solsta != mosek.solsta.integer_optimal):
                # See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
                raise Exception(f"Unexpected solution status: {solsta}")

            # Display the solution summary for quick inspection of results.
            #task.solutionsummary(mosek.streamtype.msg)

            return task.getxxslice(mosek.soltype.itg, voff_x + 0, voff_x + n)

if __name__ == '__main__':

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

    for K in range(1, n+1):
        xx = markowitz_with_card(n, k, x0, w, gamma, mu, GT, K)
        expret = sum([xx[i]*mu[i] for i in range(n)])
        print("Bound {0}:  Return = {1:.4f} x = {2}".format(K, expret, xx))

portfolio_6_factor.py

Listing 17.28 portfolio_6_factor.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      portfolio_6_factor.py
#
#  Purpose :   Implements a basic portfolio optimization model
#              of factor type.
##
import mosek
import sys
import numpy as np

if __name__ == '__main__':

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

    inf = 0.0 # This value has no significance

    # Variable offsets
    numvar = n
    voff_x = 0

    # Constraints offsets
    numcon = 1
    coff_bud = 0

    with mosek.Env() as env:
        with env.Task(0, 0) as task:
            task.set_Stream(mosek.streamtype.log, sys.stdout.write)

            # Holding variable x of length n
            # No other auxiliary variables are needed in this formulation
            task.appendvars(numvar)

            # Optionally we can give the variables names
            for j in range(0, n):
                task.putvarname(voff_x + j, "x[%d]" % (1 + j))

            # No short-selling in this model, all of x >= 0
            task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)

            # One linear constraint: total budget
            task.appendcons(1)
            task.putconname(coff_bud, "budget")
            task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n)        # e^T x
            rtemp = w + sum(x0)
            task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp)    # equals w + sum(x0)

            # Input (gamma, G_factor_T x, diag(sqrt(theta))*x) in the AFE (affine expression) storage
            # We need k+n+1 rows and we fill them in in three parts
            task.appendafes(k + n + 1)
            # 1. The first affine expression = gamma, will be specified later
            # 2. The next k expressions comprise G_factor_T*x, we add them row by row
            #    transposing the matrix G_factor on the fly 
            for i in range(0, k):
                task.putafefrow(i + 1, range(voff_x, voff_x + n), np.array(G_factor[:,i]))
            # 3. The remaining n rows contain sqrt(theta) on the diagonal
            task.putafefentrylist(range(k + 1, k + 1 + n), range(voff_x, voff_x + n), np.sqrt(theta))

            # Input the affine conic constraint (gamma, G_factor_T x, diag(sqrt(theta))*x) \in QCone
            # Add the quadratic domain of dimension k+n+1
            qdom = task.appendquadraticconedomain(k + n + 1)
            # Add the constraint
            task.appendaccseq(qdom, 0, None)    
            task.putaccname(0, "risk")

            # Objective: maximize expected return mu^T x
            task.putclist(range(voff_x, voff_x + n), mu)
            task.putobjsense(mosek.objsense.maximize)

            for gamma in gammas:
                # Specify gamma in ACC
                task.putafeg(0, gamma)

                # Dump the problem to a human readable PTF file.
                task.writedata("dump.ptf")

                # Solve the problem
                task.optimize()

                # Display solution summary for quick inspection of results.
                # In this simplified example we skip checks for problem and solution status
                task.solutionsummary(mosek.streamtype.msg)

                # Check if the interior point solution is an optimal point
                solsta = task.getsolsta(mosek.soltype.itr)
                if (solsta != mosek.solsta.optimal):
                    # See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
                    raise Exception(f"Unexpected solution status: {solsta}")

                # Retrieve results
                xx     = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
                expret = task.getprimalobj(mosek.soltype.itr)

                print(f'Expected return: {expret:.10e} Std. deviation: {gamma:.4e}')
                np.set_printoptions(precision=4)
                print(f'Optimal portfolio: {np.array(xx)}')

pow1.py

Listing 17.29 pow1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      pow1.py
#
#   Purpose: Demonstrates how to solve the problem
#
#     maximize x^0.2*y^0.8 + z^0.4 - x
#           st x + y + 0.5z = 2
#              x,y,z >= 0
##
import sys
import mosek

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()

def main():

    # Only a symbolic constant
    inf = 0.0

    # Create a task
    with mosek.Task() as task:
        # Attach a printer to the task
        task.set_Stream(mosek.streamtype.log, streamprinter)

        csub = [3, 4, 0]
        cval = [1.0, 1.0, -1.0]
        asub = [0, 1, 2]
        aval = [1.0, 1.0, 0.5]
        numvar, numcon = 5, 1         # x,y,z and 2 auxiliary variables for conic constraints

        # Append 'numcon' empty constraints.
        # The constraints will initially have no bounds.
        task.appendcons(numcon)

        # Append 'numvar' variables.
        # The variables will initially be fixed at zero (x=0).
        task.appendvars(numvar)

        # Set up the linear part of the problem
        task.putclist(csub, cval)
        task.putarow(0, asub, aval)
        task.putvarboundslice(0, numvar, [mosek.boundkey.fr] * numvar, [inf] * numvar, [inf] * numvar)
        task.putconbound(0, mosek.boundkey.fx, 2.0, 2.0)

        # Input the conic constraints
        # Append two power cone domains
        pc1 = task.appendprimalpowerconedomain(3, [0.2, 0.8])
        pc2 = task.appendprimalpowerconedomain(3, [4.0, 6.0])

        # Create data structures F,g so that
        #
        #   F * x + g = (x(0), x(1), x(3), x(2), 1.0, x(4)) 
        #
        task.appendafes(6)
        task.putafefentrylist([0, 1, 2, 3, 5],         # Rows
                              [0, 1, 3, 2, 4],         #Columns 
                              [1.0, 1.0, 1.0, 1.0, 1.0])
        task.putafeg(4, 1.0)

        # Append the two conic constraints 
        task.appendacc(pc1,                     # Domain
                       [0, 1, 2],               # Rows from F 
                       None)                    # Unused
        task.appendacc(pc2,                     # Domain
                       [3, 4, 5],               # Rows from F
                       None)                    # Unused

        # Input the objective sense (minimize/maximize)
        task.putobjsense(mosek.objsense.maximize)

        # Optimize the task
        task.optimize()

        # Print a summary containing information
        # about the solution for debugging purposes
        task.solutionsummary(mosek.streamtype.msg)
        prosta = task.getprosta(mosek.soltype.itr)
        solsta = task.getsolsta(mosek.soltype.itr)

        # Output a solution
        xx = task.getxx(mosek.soltype.itr)

        if solsta == mosek.solsta.optimal:
            print("Optimal solution: %s" % xx[0:3])
        elif solsta == mosek.solsta.dual_infeas_cer:
            print("Primal or dual infeasibility.\n")
        elif solsta == mosek.solsta.prim_infeas_cer:
            print("Primal or dual infeasibility.\n")
        elif mosek.solsta.unknown:
            print("Unknown solution status")
        else:
            print("Other solution status")


# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.code))
    if msg is not None:
        print("\t%s" % e.msg)
        sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

qcqo1.py

Listing 17.30 qcqo1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      qcqo1.py
#
#  Purpose :   Demonstrates how to solve small linear
#              optimization problem using the MOSEK Python API.
##
import sys
import mosek

# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Make a MOSEK environment
    with mosek.Env() as env:
        # Attach a printer to the environment
        env.set_Stream(mosek.streamtype.log, streamprinter)

        # Create a task
        with env.Task(0, 0) as task:
            # Attach a printer to the task
            task.set_Stream(mosek.streamtype.log, streamprinter)

            # Set up and input bounds and linear coefficients
            bkc = [mosek.boundkey.lo]
            blc = [1.0]
            buc = [inf]

            bkx = [mosek.boundkey.lo,
                   mosek.boundkey.lo,
                   mosek.boundkey.lo]
            blx = [0.0, 0.0, 0.0]
            bux = [inf, inf, inf]

            c = [0.0, -1.0, 0.0]

            asub = [[0], [0], [0]]
            aval = [[1.0], [1.0], [1.0]]

            numvar = len(bkx)
            numcon = len(bkc)
            NUMANZ = 3
            # Append 'numcon' empty constraints.
            # The constraints will initially have no bounds.
            task.appendcons(numcon)

            #Append 'numvar' variables.
            # The variables will initially be fixed at zero (x=0).
            task.appendvars(numvar)

            #Optionally add a constant term to the objective.
            task.putcfix(0.0)

            for j in range(numvar):
                # Set the linear term c_j in the objective.
                task.putcj(j, c[j])
                # Set the bounds on variable j
                # blx[j] <= x_j <= bux[j]
                task.putvarbound(j, bkx[j], blx[j], bux[j])
                # Input column j of A
                task.putacol(j,                  # Variable (column) index.
                             # Row index of non-zeros in column j.
                             asub[j],
                             aval[j])            # Non-zero Values of column j.

            for i in range(numcon):
                task.putconbound(i, bkc[i], blc[i], buc[i])

            # Set up and input quadratic objective

            qsubi = [0, 1, 2, 2]
            qsubj = [0, 1, 0, 2]
            qval = [2.0, 0.2, -1.0, 2.0]

            task.putqobj(qsubi, qsubj, qval)

            # The lower triangular part of the Q^0
            # matrix in the first constraint is specified.
            # This corresponds to adding the term
            # - x0^2 - x1^2 - 0.1 x2^2 + 0.2 x0 x2

            qsubi = [0, 1, 2, 2]
            qsubj = [0, 1, 2, 0]
            qval = [-2.0, -2.0, -0.2, 0.2]

            # put Q^0 in constraint with index 0.

            task.putqconk(0, qsubi, qsubj, qval)

            # Input the objective sense (minimize/maximize)
            task.putobjsense(mosek.objsense.minimize)

            # Optimize the task
            task.optimize()

            # Print a summary containing information
            # about the solution for debugging purposes
            task.solutionsummary(mosek.streamtype.msg)

            prosta = task.getprosta(mosek.soltype.itr)
            solsta = task.getsolsta(mosek.soltype.itr)

            # Output a solution
            xx = task.getxx(mosek.soltype.itr)

            if solsta == mosek.solsta.optimal:
                print("Optimal solution: %s" % xx)
            elif solsta == mosek.solsta.dual_infeas_cer:
                print("Primal or dual infeasibility.\n")
            elif solsta == mosek.solsta.prim_infeas_cer:
                print("Primal or dual infeasibility.\n")
            elif mosek.solsta.unknown:
                print("Unknown solution status")
            else:
                print("Other solution status")


# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    print("\t%s" % e.msg)
    sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

qo1.py

Listing 17.31 qo1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      qo1.py
#
#  Purpose :   Demonstrate how to solve a quadratic
#              optimization problem using the MOSEK Python API.
##

import sys, os, mosek

# Since the actual value of Infinity is ignored, we define it solely
# for symbolic purposes:
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Open MOSEK and create an environment and task
    # Make a MOSEK environment
    with mosek.Env() as env:
        # Attach a printer to the environment
        env.set_Stream(mosek.streamtype.log, streamprinter)
        # Create a task
        with env.Task() as task:
            task.set_Stream(mosek.streamtype.log, streamprinter)
            # Set up and input bounds and linear coefficients
            bkc = [mosek.boundkey.lo]
            blc = [1.0]
            buc = [inf]
            numvar = 3
            bkx = [mosek.boundkey.lo] * numvar
            blx = [0.0] * numvar
            bux = [inf] * numvar
            c = [0.0, -1.0, 0.0]
            asub = [[0], [0], [0]]
            aval = [[1.0], [1.0], [1.0]]

            numvar = len(bkx)
            numcon = len(bkc)

            # Append 'numcon' empty constraints.
            # The constraints will initially have no bounds.
            task.appendcons(numcon)

            # Append 'numvar' variables.
            # The variables will initially be fixed at zero (x=0).
            task.appendvars(numvar)

            for j in range(numvar):
                # Set the linear term c_j in the objective.
                task.putcj(j, c[j])
                # Set the bounds on variable j
                # blx[j] <= x_j <= bux[j]
                task.putvarbound(j, bkx[j], blx[j], bux[j])
                # Input column j of A
                task.putacol(j,                  # Variable (column) index.
                             # Row index of non-zeros in column j.
                             asub[j],
                             aval[j])            # Non-zero Values of column j.
            for i in range(numcon):
                task.putconbound(i, bkc[i], blc[i], buc[i])

            # Set up and input quadratic objective
            qsubi = [0, 1, 2, 2]
            qsubj = [0, 1, 0, 2]
            qval = [2.0, 0.2, -1.0, 2.0]

            task.putqobj(qsubi, qsubj, qval)

            # Input the objective sense (minimize/maximize)
            task.putobjsense(mosek.objsense.minimize)

            # Optimize
            task.optimize()
            # Print a summary containing information
            # about the solution for debugging purposes
            task.solutionsummary(mosek.streamtype.msg)

            prosta = task.getprosta(mosek.soltype.itr)
            solsta = task.getsolsta(mosek.soltype.itr)

            # Output a solution
            xx = task.getxx(mosek.soltype.itr)

            if solsta == mosek.solsta.optimal:
                print("Optimal solution: %s" % xx)
            elif solsta == mosek.solsta.dual_infeas_cer:
                print("Primal or dual infeasibility.\n")
            elif solsta == mosek.solsta.prim_infeas_cer:
                print("Primal or dual infeasibility.\n")
            elif mosek.solsta.unknown:
                print("Unknown solution status")
            else:
                print("Other solution status")


# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    if e.msg is not None:
        import traceback
        traceback.print_exc()
        print("\t%s" % e.msg)
    sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

reoptimization.py

Listing 17.32 reoptimization.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      reoptimization.py
#
#  Purpose:   Demonstrates how to solve a  linear
#             optimization problem using the MOSEK API
#             and modify and re-optimize the problem.
import sys
import mosek

# Since the actual value of Infinity is ignored, we define it solely
# for symbolic purposes:
inf = 0.0

def main():
    # Create a MOSEK environment
    with mosek.Env() as env:
        # Create a task
        with env.Task(0, 0) as task:
            # Bound keys for constraints
            bkc = [mosek.boundkey.up,
                   mosek.boundkey.up,
                   mosek.boundkey.up]
            # Bound values for constraints
            blc = [-inf, -inf, -inf]
            buc = [100000.0, 50000.0, 60000.0]
            # Bound keys for variables
            bkx = [mosek.boundkey.lo,
                   mosek.boundkey.lo,
                   mosek.boundkey.lo]
            # Bound values for variables
            blx = [0.0, 0.0, 0.0]
            bux = [+inf, +inf, +inf]
            # Objective coefficients
            csub = [0, 1, 2]
            cval = [1.5, 2.5, 3.0]
            # We input the A matrix column-wise
            # asub contains row indexes
            asub = [0, 1, 2,
                    0, 1, 2,
                    0, 1, 2]
            # acof contains coefficients
            acof = [2.0, 3.0, 2.0,
                    4.0, 2.0, 3.0,
                    3.0, 3.0, 2.0]
            # aptrb and aptre contains the offsets into asub and acof where
            # columns start and end respectively
            aptrb = [0, 3, 6]
            aptre = [3, 6, 9]

            numvar = len(bkx)
            numcon = len(bkc)

            # Append the constraints
            task.appendcons(numcon)

            # Append the variables.
            task.appendvars(numvar)

            # Input objective
            task.putcfix(0.0)
            task.putclist(csub, cval)

            # Put constraint bounds
            task.putconboundslice(0, numcon, bkc, blc, buc)

            # Put variable bounds
            task.putvarboundslice(0, numvar, bkx, blx, bux)

            # Input A non-zeros by columns
            for j in range(numvar):
                ptrb, ptre = aptrb[j], aptre[j]
                task.putacol(j,
                             asub[ptrb:ptre],
                             acof[ptrb:ptre])

            # Input the objective sense (minimize/maximize)
            task.putobjsense(mosek.objsense.maximize)

            # Optimize the task
            task.optimize()

            # Output a solution
            xx = task.getsolutionslice(mosek.soltype.bas,
                                       mosek.solitem.xx,
                                       0, numvar)
            print("xx = {}".format(xx))

            ################# Make a change to the A matrix #############
            task.putaij(0, 0, 3.0)
            task.optimize()
            # Output a solution
            xx = task.getsolutionslice(mosek.soltype.bas,
                                       mosek.solitem.xx,
                                       0, numvar)
            print("xx = {}".format(xx))

            ################### Add a new variable ######################
            task.appendvars(1)
            numvar+=1

            # Set bounds on new varaible
            task.putvarbound(task.getnumvar() - 1,
                             mosek.boundkey.lo,
                             0,
                             +inf)

            # Change objective
            task.putcj(task.getnumvar() - 1, 1.0)

            # Put new values in the A matrix
            acolsub = [0, 2]
            acolval = [4.0, 1.0]

            task.putacol(task.getnumvar() - 1, # column index
                         acolsub,
                         acolval)
            # Change optimizer to simplex free and reoptimize
            task.putintparam(mosek.iparam.optimizer,
                             mosek.optimizertype.free_simplex)
            task.optimize()
            # Output a solution
            xx = task.getsolutionslice(mosek.soltype.bas,
                                       mosek.solitem.xx,
                                       0, numvar)
            print("xx = {}".format(xx))

            ############# Add a new constraint #######################
            task.appendcons(1)
            numcon+=1

            # Set bounds on new constraint
            task.putconbound(task.getnumcon() - 1,
                             mosek.boundkey.up, -inf, 30000)

            # Put new values in the A matrix
            arowsub = [0, 1, 2, 3]
            arowval = [1.0, 2.0, 1.0, 1.0]

            task.putarow(task.getnumcon() - 1, # row index
                         arowsub,
                         arowval)

            task.optimize()
            # Output a solution
            xx = task.getsolutionslice(mosek.soltype.bas,
                                       mosek.solitem.xx,
                                       0, numvar)
            print("xx = {}".format(xx))

            ############# Change constraint bounds #######################
            newbkc = [mosek.boundkey.up] * numcon
            newblc = [-inf] * numcon
            newbuc = [ 80000, 40000, 50000, 22000 ]

            task.putconboundslice(0, numcon, newbkc, newblc, newbuc)

            task.optimize()

            # Output a solution
            xx = task.getsolutionslice(mosek.soltype.bas,
                                       mosek.solitem.xx,
                                       0, numvar)
            print("xx = {}".format(xx))



# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    print("\t%s" % e.msg)
    sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

response.py

Listing 17.33 response.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      response.py
#
#  Purpose :   This example demonstrates proper response handling
#              for problems solved with the interior-point optimizers.
#
##
import mosek
import sys

# A log message
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()

def main(args):
  filename = args[0] if len(args) >= 1 else "../data/cqo1.mps"

  try:
    # Create environment and task
    with mosek.Env() as env:
      with env.Task(0, 0) as task:
        # (Optional) set a log stream
        # task.set_Stream(mosek.streamtype.log, streamprinter)

        # (Optional) uncomment to see what happens when solution status is unknown
        # task.putintparam(mosek.iparam.intpnt_max_iterations, 1)

        # In this example we read data from a file
        task.readdata(filename)

        # Optimize
        trmcode = task.optimize()
        task.solutionsummary(mosek.streamtype.log)

        # We expect solution status OPTIMAL
        solsta = task.getsolsta(mosek.soltype.itr)

        if solsta == mosek.solsta.optimal:
          # Optimal solution. Fetch and print it.
          print("An optimal interior-point solution is located.")
          numvar = task.getnumvar()
          xx = task.getxx(mosek.soltype.itr)
          for i in range(numvar): 
            print("x[{0}] = {1}".format(i, xx[i]))

        elif solsta == mosek.solsta.dual_infeas_cer:
          print("Dual infeasibility certificate found.")

        elif solsta == mosek.solsta.prim_infeas_cer:
          print("Primal infeasibility certificate found.")
        
        elif solsta == mosek.solsta.unknown:
          # The solutions status is unknown. The termination code
          # indicates why the optimizer terminated prematurely.
          print("The solution status is unknown.")
          symname, desc = mosek.Env.getcodedesc(trmcode)
          print("   Termination code: {0} {1}".format(symname, desc))

        else:
          print("An unexpected solution status {0} is obtained.".format(str(solsta)))

  except mosek.Error as e:
      print("Unexpected error ({0}) {1}".format(e.errno, e.msg))

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

sdo1.py

Listing 17.34 sdo1.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      sdo1.py
#
#  Purpose :   Demonstrates how to solve a small mixed semidefinite and conic quadratic
#              optimization problem using the MOSEK Python API.
##
import sys
import mosek

# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Create a task object and attach log stream printer
    with mosek.Task() as task:
        task.set_Stream(mosek.streamtype.log, streamprinter)

        # Bound keys for constraints
        bkc = [mosek.boundkey.fx,
               mosek.boundkey.fx]

        # Bound values for constraints
        blc = [1.0, 0.5]
        buc = [1.0, 0.5]

        # Below is the sparse representation of the A
        # matrix stored by row.
        asub = [[0],
                [1, 2]]
        aval = [[1.0],
                [1.0, 1.0]]

        barci = [0, 1, 1, 2, 2]
        barcj = [0, 0, 1, 1, 2]
        barcval = [2.0, 1.0, 2.0, 1.0, 2.0]

        barai = [[0, 1, 2],
                 [0, 1, 2, 1, 2, 2]]
        baraj = [[0, 1, 2],
                 [0, 0, 0, 1, 1, 2]]
        baraval = [[1.0, 1.0, 1.0],
                   [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]]

        numvar = 3
        numcon = len(bkc)
        BARVARDIM = [3]

        # Append 'numvar' variables.
        # The variables will initially be fixed at zero (x=0).
        task.appendvars(numvar)

        # Append 'numcon' empty constraints.
        # The constraints will initially have no bounds.
        task.appendcons(numcon)

        # Append matrix variables of sizes in 'BARVARDIM'.
        # The variables will initially be fixed at zero.
        task.appendbarvars(BARVARDIM)

        # Set the linear term c_0 in the objective.
        task.putcj(0, 1.0)

        for j in range(numvar):
            # Set the bounds on variable j
            # blx[j] <= x_j <= bux[j]
            task.putvarbound(j, mosek.boundkey.fr, -inf, +inf)

        for i in range(numcon):
            # Set the bounds on constraints.
            # blc[i] <= constraint_i <= buc[i]
            task.putconbound(i, bkc[i], blc[i], buc[i])

            # Input row i of A
            task.putarow(i,                  # Constraint (row) index.
                         asub[i],            # Column index of non-zeros in constraint i.
                         aval[i])            # Non-zero values of row i.

        # Add the quadratic cone constraint
        task.appendafes(3)
        # Diagonal F matrix
        task.putafefentrylist(range(3), range(3), [1.0]*3)
        task.appendaccseq(task.appendquadraticconedomain(3), 0, None)
        
        symc = task.appendsparsesymmat(BARVARDIM[0],
                                       barci,
                                       barcj,
                                       barcval)

        syma0 = task.appendsparsesymmat(BARVARDIM[0],
                                        barai[0],
                                        baraj[0],
                                        baraval[0])

        syma1 = task.appendsparsesymmat(BARVARDIM[0],
                                        barai[1],
                                        baraj[1],
                                        baraval[1])

        task.putbarcj(0, [symc], [1.0])

        task.putbaraij(0, 0, [syma0], [1.0])
        task.putbaraij(1, 0, [syma1], [1.0])

        # Input the objective sense (minimize/maximize)
        task.putobjsense(mosek.objsense.minimize)

        # Solve the problem and print summary
        task.optimize()
        task.solutionsummary(mosek.streamtype.msg)

        # Get status information about the solution
        prosta = task.getprosta(mosek.soltype.itr)
        solsta = task.getsolsta(mosek.soltype.itr)

        if (solsta == mosek.solsta.optimal):
            xx = task.getxx(mosek.soltype.itr)
            barx = task.getbarxj(mosek.soltype.itr, 0)

            print("Optimal solution:\nx=%s\nbarx=%s" % (xx, barx))
        elif (solsta == mosek.solsta.dual_infeas_cer or
              solsta == mosek.solsta.prim_infeas_cer):
            print("Primal or dual infeasibility certificate found.\n")
        elif solsta == mosek.solsta.unknown:
            print("Unknown solution status")
        else:
            print("Other solution status")


# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    if e.msg is not None:
        print("\t%s" % e.msg)
        sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

sdo2.py

Listing 17.35 sdo2.py Click here to download.
#
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      sdo2.py
#
#  Purpose :   Solves the semidefinite problem with two symmetric variables:
#
#                 min   <C1,X1> + <C2,X2>
#                 st.   <A1,X1> + <A2,X2> = b
#                             (X2)_{1,2} <= k
#                
#                 where X1, X2 are symmetric positive semidefinite,
#
#                 C1, C2, A1, A2 are assumed to be constant symmetric matrices,
#                 and b, k are constants.

from mosek import *
import sys

# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()
                
# Sample data in sparse lower-triangular triplet form
C1_k = [0, 2]
C1_l = [0, 2]
C1_v = [1, 6]
A1_k = [0, 2, 2]
A1_l = [0, 0, 2]
A1_v = [1, 1, 2]
C2_k = [0, 1, 1, 2]
C2_l = [0, 0, 1, 2]
C2_v = [1, -3, 2, 1]
A2_k = [1, 1, 3]
A2_l = [0, 1, 3]
A2_v = [1, -1, -3]
b = 23
k = -3

# Make mosek environment
with Env() as env:

    # Create a task object and attach log stream printer
    with env.Task(0, 0) as task:
        # Set log handler for debugging ootput
        task.set_Stream(streamtype.log, streamprinter)

        # Append two symmetric variables of dimension 3, 4
        barvardims = [3, 4]
        task.appendbarvars(barvardims)

        # Semidefinite part of objective function
        task.putbarcblocktriplet(
            [0]*len(C1_v) + [1]*len(C2_v), # Which SDP variable (j)
            C1_k + C2_k,                   # Entries: (k,l)->v
            C1_l + C2_l,
            C1_v + C2_v,
            )

        # Append two constraints
        task.appendcons(2)

        # First constraint (equality)
        task.putbarablocktriplet(
            [0]*(len(A1_v)+len(A2_v)),     # Which constraint (i = 0)
            [0]*len(A1_v) + [1]*len(A2_v), # Which SDP variable (j)
            A1_k + A2_k,                   # Entries: (k,l)->v
            A1_l + A2_l,
            A1_v + A2_v,
            )

        # Second constraint (X2)_{1,2} <= k
        task.putbarablocktriplet(
            [1],                           # Which constraint (i = 1)
            [1],                           # Which SDP variable (j = 1)
            [1], [0], [0.5]                # Entries: (k,l)->v
            )

        # Set bounds for constraints
        task.putconboundlist([0,1], [boundkey.fx, boundkey.up],
                                    [b, -inf],
                                    [b, k])

        # Write the problem for human inspection
        task.writedata("test.ptf")

        # Optimize
        task.optimize()
        task.solutionsummary(streamtype.msg)

        # Get status information about the solution
        solsta = task.getsolsta(soltype.itr)

        if solsta == solsta.optimal:
            # Assuming the optimization succeeded read solution
            print("Solution (lower-triangular part vectorized): ")
            for i in range(2):
                X = task.getbarxj(soltype.itr, i)
                print("X{i} = {X}".format(i=i, X=X))

        elif (solsta == solsta.dual_infeas_cer or
              solsta == solsta.prim_infeas_cer):
            print("Primal or dual infeasibility certificate found.\n")
        elif solsta == solsta.unknown:
            print("Unknown solution status")
        else:
            print("Other solution status")

sdo_lmi.py

Listing 17.36 sdo_lmi.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      sdo_lmi.py
#
#  Purpose :   To solve a problem with an LMI and an affine conic constrained problem with a PSD term
#
#                 minimize    Tr [1, 0; 0, 1]*X + x(1) + x(2) + 1
#
#                 subject to  Tr [0, 1; 1, 0]*X - x(1) - x(2) >= 0
#                             x(1) [0, 1; 1, 3] + x(2) [3, 1; 1, 0] - [1, 0; 0, 1] >> 0
#                             X >> 0
import sys
from numpy import sqrt
import mosek

# Since the value of infinity is ignored, we define it solely
# for symbolic purposes
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Create a task object and attach log stream printer
    with mosek.Task() as task:
        task.set_Stream(mosek.streamtype.log, streamprinter)

        # Below is the sparse triplet representation of the F matrix.
        afeidx = [0, 0, 1, 2, 2, 3]
        varidx = [0, 1, 1, 0, 1, 0]
        f_val  = [-1, -1, 3, sqrt(2), sqrt(2), 3]
        g      = [0, -1, 0, -1]

        barcj = [0, 0]
        barck = [0, 1]
        barcl = [0, 1]
        barcv = [1, 1]

        barfi = [0,0]
        barfj = [0,0]
        barfk = [0,1]
        barfl = [0,0]
        barfv = [0,1]

        numvar = 2
        numafe = 4
        BARVARDIM = [2]

        # Append 'numvar' variables.
        # The variables will initially be fixed at zero (x=0).
        task.appendvars(numvar)

        # Append 'numafe' empty affine expressions.
        task.appendafes(numafe)

        # Append matrix variables of sizes in 'BARVARDIM'.
        # The variables will initially be fixed at zero.
        task.appendbarvars(BARVARDIM)

        # Set the linear terms in the objective.
        task.putcj(0, 1.0)
        task.putcj(1, 1.0)
        task.putcfix(1.0)
        task.putbarcblocktriplet(barcj, barck, barcl, barcv)

        for j in range(numvar):
            # Set the bounds on variable j
            # blx[j] <= x_j <= bux[j]
            task.putvarbound(j, mosek.boundkey.fr, -inf, +inf)

        # Set up the F matrix of the problem
        task.putafefentrylist(afeidx, varidx, f_val)
        # Set up the g vector of the problem
        task.putafegslice(0, numafe, g)
        task.putafebarfblocktriplet(barfi, barfj, barfk, barfl, barfv)

        # Append R+ domain and the corresponding ACC
        task.appendacc(task.appendrplusdomain(1), [0], None)
        # Append SVEC_PSD domain and the corresponding ACC
        task.appendacc(task.appendsvecpsdconedomain(3), [1,2,3], None)

        # Input the objective sense (minimize/maximize)
        task.putobjsense(mosek.objsense.minimize)

        # Solve the problem and print summary
        task.optimize()
        task.solutionsummary(mosek.streamtype.msg)

        # Get status information about the solution
        prosta = task.getprosta(mosek.soltype.itr)
        solsta = task.getsolsta(mosek.soltype.itr)

        if (solsta == mosek.solsta.optimal):
            xx = task.getxx(mosek.soltype.itr)
            barx = task.getbarxj(mosek.soltype.itr, 0)

            print("Optimal solution:\nx=%s\nbarx=%s" % (xx, barx))
        elif (solsta == mosek.solsta.dual_infeas_cer or
              solsta == mosek.solsta.prim_infeas_cer):
            print("Primal or dual infeasibility certificate found.\n")
        elif solsta == mosek.solsta.unknown:
            print("Unknown solution status")
        else:
            print("Other solution status")


# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    if e.msg is not None:
        print("\t%s" % e.msg)
        sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

sensitivity.py

Listing 17.37 sensitivity.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      sensitivity.py
#
#  Purpose :   To demonstrate how to perform sensitivity
#              analysis from the API on a small problem:
#
#              minimize
#
#              obj: +1 x11 + 2 x12 + 5 x23 + 2 x24 + 1 x31 + 2 x33 + 1 x34
#              st
#              c1:   +  x11 +   x12                                           <= 400
#              c2:                  +   x23 +   x24                           <= 1200
#              c3:                                  +   x31 +   x33 +   x34   <= 1000
#              c4:   +  x11                         +   x31                   = 800
#              c5:          +   x12                                           = 100
#              c6:                  +   x23                 +   x33           = 500
#              c7:                          +   x24                 +   x34   = 500
#
#              The example uses basis type sensitivity analysis.
##

import sys
import mosek

# Since the actual value of Infinity is ignores, we define it solely
# for symbolic purposes:
inf = 0.0

# Define a stream printer to grab output from MOSEK
def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()


def main():
    # Create a MOSEK environment
    with mosek.Env() as env:
        # Attach a printer to the environment
        env.set_Stream(mosek.streamtype.log, streamprinter)

        # Create a task
        with env.Task(0, 0) as task:
            # Attach a printer to the task
            task.set_Stream(mosek.streamtype.log, streamprinter)

            # Set up data
            bkc = [mosek.boundkey.up, mosek.boundkey.up,
                   mosek.boundkey.up, mosek.boundkey.fx,
                   mosek.boundkey.fx, mosek.boundkey.fx,
                   mosek.boundkey.fx]
            blc = [-inf, -inf, -inf, 800., 100., 500., 500.]
            buc = [400., 1200., 1000., 800., 100., 500., 500.]

            bkx = [mosek.boundkey.lo, mosek.boundkey.lo,
                   mosek.boundkey.lo, mosek.boundkey.lo,
                   mosek.boundkey.lo, mosek.boundkey.lo,
                   mosek.boundkey.lo]
            c = [1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 1.0]
            blx = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
            bux = [inf, inf, inf, inf, inf, inf, inf]

            ptrb = [0, 2, 4, 6, 8, 10, 12]
            ptre = [2, 4, 6, 8, 10, 12, 14]
            sub = [0, 3, 0, 4, 1, 5, 1, 6, 2, 3, 2, 5, 2, 6]

            val = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
                   1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

            numcon = len(bkc)
            numvar = len(bkx)
            numanz = len(val)

            # Input linear data
            task.inputdata(numcon, numvar,
                           c, 0.0,
                           ptrb, ptre, sub, val,
                           bkc, blc, buc,
                           bkx, blx, bux)
            # Set objective sense
            task.putobjsense(mosek.objsense.minimize)

            # Optimize
            task.optimize()

            # Analyze upper bound on c1 and the equality constraint on c4
            subi = [0, 3]
            marki = [mosek.mark.up, mosek.mark.up]

            # Analyze lower bound on the variables x12 and x31
            subj = [1, 4]
            markj = [mosek.mark.lo, mosek.mark.lo]

            (leftpricei, rightpricei, leftrangei, rightrangei, 
             leftpricej, rightpricej, leftrangej, rightrangej) = task.primalsensitivity(subi,
                                                                                        marki,
                                                                                        subj,
                                                                                        markj)

            print('Results from sensitivity analysis on bounds:')
            print('\tleftprice  | rightprice | leftrange  | rightrange ')
            print('For constraints:')

            for i in range(2):
                print('\t%10f   %10f   %10f   %10f' % (leftpricei[i],
                                                       rightpricei[i],
                                                       leftrangei[i],
                                                       rightrangei[i]))

            print('For variables:')
            for i in range(2):
                print('\t%10f   %10f   %10f   %10f' % (leftpricej[i],
                                                       rightpricej[i],
                                                       leftrangej[i],
                                                       rightrangej[i]))

            subc = [2, 5]

            leftprice, rightprice, leftrange, rightrange = task.dualsensitivity(subc)

            print('Results from sensitivity analysis on objective coefficients:')

            for i in range(2):
                print('\t%10f   %10f   %10f   %10f' % (leftprice[i],
                                                       rightprice[i],
                                                       leftrange[i],
                                                       rightrange[i]))

    return None

# call the main function
try:
    main()
except mosek.MosekException as e:
    print("ERROR: %s" % str(e.errno))
    if e.msg is not None:
        print("\t%s" % e.msg)
    sys.exit(1)
except:
    import traceback
    traceback.print_exc()
    sys.exit(1)

simple.py

Listing 17.38 simple.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      simple.py
#
#  Purpose :   Demonstrates a very simple example using MOSEK by
#              reading a problem file, solving the problem and
#              writing the solution to a file.
##
import mosek
import sys

def streamprinter(msg):
    sys.stdout.write(msg)
    sys.stdout.flush()

if len(sys.argv) <= 1:
    print("Missing argument, syntax is:")
    print("  simple inputfile [ solutionfile ]")
else:
    # Create the mosek environment.
    with mosek.Env() as env:

        # Create a task object linked with the environment env.
        # We create it with 0 variables and 0 constraints initially,
        # since we do not know the size of the problem.
        with env.Task(0, 0) as task:
            task.set_Stream(mosek.streamtype.log, streamprinter)

            # We assume that a problem file was given as the first command
            # line argument (received in `argv')

            task.readdata(sys.argv[1])

            # Solve the problem
            task.optimize()

            # Print a summary of the solution
            task.solutionsummary(mosek.streamtype.log)

            # If an output file was specified, save problem to a file
            if len(sys.argv) >= 3:
                # If using OPF format, these parameters will specify what to include in output
                task.putintparam(mosek.iparam.opf_write_solutions, mosek.onoffkey.on)
                task.putintparam(mosek.iparam.opf_write_problem, mosek.onoffkey.on)
                task.putintparam(mosek.iparam.opf_write_hints, mosek.onoffkey.off)
                task.putintparam(mosek.iparam.opf_write_parameters, mosek.onoffkey.off)

                task.writedata(sys.argv[2])

solutionquality.py

Listing 17.39 solutionquality.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      solutionquality.py
#
#  Purpose :   To demonstrate how to examine the quality of a solution.
##
import sys
import mosek

def streamprinter(msg):
    sys.stdout.write(msg)
    sys.stdout.flush()

if len(sys.argv) <= 1:
    print("Missing argument, syntax is:")
    print("  solutionquality inputfile")
else:
    try:
        # Create the mosek environment.
        with mosek.Env() as env:
            # Create a task object linked with the environment env.
            # We create it with 0 variables and 0 constraints initially,
            # since we do not know the size of the problem.
            with env.Task(0, 0) as task:

                task.set_Stream(mosek.streamtype.log, streamprinter)

                # We assume that a problem file was given as the first command
                # line argument (received in `argv')
                task.readdata(sys.argv[1])

                # Solve the problem
                task.optimize()

                # Print a summary of the solution
                task.solutionsummary(mosek.streamtype.log)

                whichsol = mosek.soltype.bas

                solsta = task.getsolsta(whichsol)

                pobj, pviolcon, pviolvar, pviolbarvar, pviolcones, pviolitg, \
                dobj, dviolcon, dviolvar, dviolbarvar, dviolcones = \
                    task.getsolutioninfo(whichsol)

                if solsta in [mosek.solsta.optimal]:

                    abs_obj_gap = abs(dobj - pobj)
                    rel_obj_gap = abs_obj_gap / \
                        (1.0 + min(abs(pobj), abs(dobj)))
                    max_primal_viol = max(pviolcon, pviolvar)
                    max_primal_viol = max(max_primal_viol, pviolbarvar)
                    max_primal_viol = max(max_primal_viol, pviolcones)

                    max_dual_viol = max(dviolcon, dviolvar)
                    max_dual_viol = max(max_dual_viol, dviolbarvar)
                    max_dual_viol = max(max_dual_viol, dviolcones)

                    # Assume the application needs the solution to be within
                    #    1e-6 ofoptimality in an absolute sense. Another approach
                    #   would be looking at the relative objective gap

                    print("\n\n")
                    print("Customized solution information.\n")
                    print("  Absolute objective gap: %e\n" % abs_obj_gap)
                    print("  Relative objective gap: %e\n" % rel_obj_gap)
                    print("  Max primal violation  : %e\n" % max_primal_viol)
                    print("  Max dual violation    : %e\n" % max_dual_viol)

                    accepted = True

                    if rel_obj_gap > 1e-6:
                        print("Warning: The relative objective gap is LARGE.")
                        accepted = False

                    # We will accept a primal infeasibility of 1e-8 and
                    # dual infeasibility of 1e-6. These number should chosen problem
                    # dependent.
                    if max_primal_viol > 1e-8:
                        print("Warning: Primal violation is too LARGE")
                        accepted = False

                    if max_dual_viol > 1e-6:
                        print("Warning: Dual violation is too LARGE.")
                        accepted = False

                    if accepted:

                        numvar = task.getnumvar()
                        print("Optimal primal solution")
                        x = task.getxx(whichsol)
                        for j in range(numvar):
                            print("x[%d]: %e\n" % (j, x[j]))

                    else:
                        #Print detailed information about the solution
                        task.analyzesolution(mosek.streamtype.log, whichsol)

                elif solsta in [mosek.solsta.dual_infeas_cer, mosek.solsta.prim_infeas_cer]:

                    print("Primal or dual infeasibility certificate found.")

                elif solsta == mosek.solsta.unkwown:
                    print("The status of the solution is unknown.")
                else:
                    print("Other solution status")

    except mosek.Error as e:
        print(e)

solvebasis.py

Listing 17.40 solvebasis.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      solvebasis.py
#
#  Purpose :   To demonstrate the usage of
#              MSK_solvewithbasis on the problem:
#
#              maximize  x0 + x1
#              st.
#                      x0 + 2.0 x1 <= 2
#                      x0  +    x1 <= 6
#                      x0 >= 0, x1>= 0
#
#               The problem has the slack variables
#               xc0, xc1 on the constraints
#               and the variabels x0 and x1.
#
#               maximize  x0 + x1
#               st.
#                  x0 + 2.0 x1 -xc1       = 2
#                  x0  +    x1       -xc2 = 6
#                  x0 >= 0, x1>= 0,
#                  xc1 <=  0 , xc2 <= 0
##
import mosek
import sys

def streamprinter(text):
    sys.stdout.write(text)
    sys.stdout.flush()

def main():
    numcon = 2
    numvar = 2

    # Since the value infinity is never used, we define
    # 'infinity' symbolic purposes only
    infinity = 0

    c = [1.0, 1.0]
    ptrb = [0, 2]
    ptre = [2, 3]
    asub = [0, 1,
            0, 1]
    aval = [1.0, 1.0,
            2.0, 1.0]
    bkc = [mosek.boundkey.up,
           mosek.boundkey.up]

    blc = [-infinity,
           -infinity]
    buc = [2.0,
           6.0]

    bkx = [mosek.boundkey.lo,
           mosek.boundkey.lo]
    blx = [0.0,
           0.0]

    bux = [+infinity,
           +infinity]
    w1 = [2.0, 6.0]
    w2 = [1.0, 0.0]

    try:
        with mosek.Env() as env:
            with env.Task(0, 0) as task:
                task.set_Stream(mosek.streamtype.log, streamprinter)
                task.inputdata(numcon, numvar,
                               c,
                               0.0,
                               ptrb,
                               ptre,
                               asub,
                               aval,
                               bkc,
                               blc,
                               buc,
                               bkx,
                               blx,
                               bux)
                task.putobjsense(mosek.objsense.maximize)
                r = task.optimize()
                if r != mosek.rescode.ok:
                    print("Mosek warning:", r)

                basis = task.initbasissolve()

                #List basis variables corresponding to columns of B
                varsub = [0, 1]

                for i in range(numcon):
                    if basis[varsub[i]] < numcon:
                        print("Basis variable no %d is xc%d" % (i, basis[i]))
                    else:
                        print("Basis variable no %d is x%d" %
                              (i, basis[i] - numcon))

                # solve Bx = w1
                # varsub contains index of non-zeros in b.
                #  On return b contains the solution x and
                # varsub the index of the non-zeros in x.
                nz = 2

                nz = task.solvewithbasis(False, nz, varsub, w1)
                print("nz = %s" % nz)
                print("Solution to Bx = w1:")

                for i in range(nz):
                    if basis[varsub[i]] < numcon:
                        print("xc %s = %s" % (basis[varsub[i]], w1[varsub[i]]))
                    else:
                        print("x%s = %s" %
                              (basis[varsub[i]] - numcon, w1[varsub[i]]))

                # Solve B^Tx = w2
                nz = 1
                varsub[0] = 0

                nz = task.solvewithbasis(True, nz, varsub, w2)

                print("Solution to B^Tx = w2:")

                for i in range(nz):
                    if basis[varsub[i]] < numcon:
                        print("xc %s = %s" % (basis[varsub[i]], w2[varsub[i]]))
                    else:
                        print("x %s = %s" %
                              (basis[varsub[i]] - numcon, w2[varsub[i]]))
    except Exception as e:
        print(e)

if __name__ == '__main__':
    main()

solvelinear.py

Listing 17.41 solvelinear.py Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      solvelinear.py
#
#  Purpose :   To demonstrate the usage of MSK_solvewithbasis
#              when solving the linear system:
#
#                1.0  x1             = b1
#               -1.0  x0  +  1.0  x1 = b2
#
#              with two different right hand sides
#
#                 b = (1.0, -2.0)
#
#              and
#
#                 b = (7.0, 0.0)
##
import mosek

def setup(task,
          aval,
          asub,
          ptrb,
          ptre,
          numvar,
          basis):
    # Since the value infinity is never used, we define
    # 'infinity' symbolic purposes only
    infinity = 0

    skx = [mosek.stakey.bas] * numvar
    skc = [mosek.stakey.fix] * numvar

    task.appendvars(numvar)
    task.appendcons(numvar)

    for i in range(len(asub)):
        task.putacol(i, asub[i], aval[i])

    for i in range(numvar):
        task.putconbound(i, mosek.boundkey.fx, 0.0, 0.0)

    for i in range(numvar):
        task.putvarbound(i,
                         mosek.boundkey.fr,
                         -infinity,
                         infinity)

    # Define a basic solution by specifying
    # status keys for variables & constraints.
    task.deletesolution(mosek.soltype.bas);

    task.putskcslice(mosek.soltype.bas, 0, numvar, skc);
    task.putskxslice(mosek.soltype.bas, 0, numvar, skx);

    task.initbasissolve(basis);


def main():
    numcon = 2
    numvar = 2

    aval = [[-1.0],
            [1.0, 1.0]]
    asub = [[1],
            [0, 1]]

    ptrb = [0, 1]
    ptre = [1, 3]

    with mosek.Task() as task:
        # Directs the log task stream to the user specified
        # method task_msg_obj.streamCB
        task.set_Stream(mosek.streamtype.log,
                        lambda msg: sys.stdout.write(msg))
        # Put A matrix and factor A.
        # Call this function only once for a given task.

        basis = [0] * numvar
        b = [0.0, -2.0]
        bsub = [0, 1]

        setup(task,
              aval,
              asub,
              ptrb,
              ptre,
              numvar,
              basis)

        # now solve rhs
        b = [1, -2]
        bsub = [0, 1]
        nz = task.solvewithbasis(False, 2, bsub, b)
        print("\nSolution to Bx = b:\n")

        # Print solution and show correspondents
        # to original variables in the problem
        for i in range(nz):
            if basis[bsub[i]] < numcon:
                print("This should never happen")
            else:
                print("x%d = %d" % (basis[bsub[i]] - numcon, b[bsub[i]]))

        b[0] = 7
        bsub[0] = 0

        nz = task.solvewithbasis(False, 1, bsub, b)

        print("\nSolution to Bx = b:\n")
        # Print solution and show correspondents
        # to original variables in the problem
        for i in range(nz):
            if basis[bsub[i]] < numcon:
                print("This should never happen")
            else:
                print("x%d = %d" % (basis[bsub[i]] - numcon, b[bsub[i]]))

if __name__ == "__main__":
    try:
        main()
    except:
        import traceback
        traceback.print_exc()

sparsecholesky.py

Listing 17.42 sparsecholesky.py Click here to download.
##
#  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File:      sparsecholesky.py
#
#  Purpose: Demonstrate the sparse Cholesky factorization.
##

from __future__ import print_function
import mosek
import numpy

# Prints out a Cholesky factor presented in sparse form
def printsparse(n, perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc):
    print("P = ", perm)
    print("diag(D) = ", numpy.array(diag))
    print("L=")
    l = [[0.0 for j in range(n)] for i in range(n)]
    for j in range(n):
        for i in range(lptrc[j], lptrc[j] + lnzc[j]):
            l[lsubc[i]][j] = lvalc[i]
    print(numpy.array(l))


# Create the mosek environment.
with mosek.Env() as env:
    # The matrix from the manual
    # Observe that anzc, aptrc, asubc and avalc only specify the lower
    # triangular part.
    n = 4
    anzc = [4, 1, 1, 1]
    asubc = [0, 1, 2, 3, 1, 2, 3]
    aptrc = [0, 4, 5, 6]
    avalc = [4.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    b = [13.0, 3.0, 4.0, 5.0]

    try:
        perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc = env.computesparsecholesky(
            0,      # Mosek chooses number of threads
            1,      # Use reordering heuristic
            1.0e-14,# Singularity tolerance
            anzc, aptrc, asubc, avalc)

        printsparse(n, perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc)

        x = [b[p] for p in perm]   # Permuted b is stored as x.

        # Compute inv(L)*x.
        env.sparsetriangularsolvedense(mosek.transpose.no,
                                       lnzc, lptrc, lsubc, lvalc, x)

        # Compute inv(L^T)*x.
        env.sparsetriangularsolvedense(mosek.transpose.yes,
                                       lnzc, lptrc, lsubc, lvalc, x)

        print("\nSolution Ax=b: x = ", numpy.array(
            [x[j] for i in range(n) for j in range(n) if perm[j] == i]), "\n")
    except:
        raise

    #Example 2 - singular A
    n = 3
    anzc = [3, 2, 1]
    asubc = [0, 1, 2, 1, 2, 2]
    aptrc = [0, 3, 5]
    avalc = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

    try:
        perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc = env.computesparsecholesky(
            0,      # Mosek chooses number of threads
            1,      # Use reordering heuristic
            1.0e-14,# Singularity tolerance
            anzc, aptrc, asubc, avalc)

        printsparse(n, perm, diag, lnzc, lptrc, lensubnval, lsubc, lvalc)
    except:
        raise