6. Advanced API tutorial


This chapter provides information about additional problem classes and functionality provided in the Python API.

6.1. Linear network flow problems

Network flow problems are a special class of linear optimization problems which has many applications. A network consists of a set of points connected by a set of lines. Usually the points and lines are called nodes and arcs. Arcs may have an direction on them. The network is directed if all arcs are directed. The class of network flow problems is defined as follows.
Let [[MathCmd 119]] be a directed network of nodes [[MathCmd 120]] and arcs [[MathCmd 121]]. Associated with every arc [[MathCmd 122]] is a cost [[MathCmd 123]] and a capacity [[MathCmd 124]]. Moreover, associated with each node [[MathCmd 125]] in the network is a lower limit [[MathCmd 126]] and an upper limit [[MathCmd 127]] on the demand (supply) of the node. The minimum cost of a network flow problem can be stated as follows:

\begin{math}\nonumber{}\begin{array}{lcccccl}\nonumber{}\mbox{minimize} &  &  & \sum \limits _{{(i,j)\in{}\mathcal{A}}}c_{{ij}}x_{{ij}} &  &  & \\\nonumber{}\mbox{subject to} & l^{c}_{{i}} & \leq{} & \sum \limits _{{\lbrace{}j:(i,j)\in{}\mathcal{A}\rbrace{}}}x_{{ij}}-\sum \limits _{{\lbrace{}j:(j,i)\in{}\mathcal{A}\rbrace{}}}x_{{ji}} & \leq{} & u^{c}_{{i}} & \forall i\in{}\mathcal{N},\\\nonumber{} & l^{x}_{{ij}} & \leq{} & x_{{ij}} & \leq{} & u^{x}_{{ij}} & \forall (i,j)\in{}\mathcal{A}.\end{array}\end{math} (6.1.1)

A classical example of a network flow problem is the transportation problem where the objective is to distribute goods from warehouses to customers at lowest possible total cost, see [7] for a detailed application reference.

The above graph formulation of the network flow problem implies the structural properties. Each variable appears in exactly two constraints with a numerical value of either [[MathCmd 129]] or [[MathCmd 130]].

It is well-known that problems with network flow structure can be solved efficiently with a specialized version of the simplex method. MOSEK includes such a network simplex implementation which can be called either directly using Task.netoptimize or indirectly by letting the standard simplex optimizer extract the embedded network. This section shows how to solve a network problem by a direct call to Task.netoptimize. For further details on how to exploit embedded network in the standard simplex optimizer, see Section 8.3.1.

6.1.1. A linear network flow problem example

The following is an example of a linear network optimization problem:

\begin{math}\nonumber{}\begin{array}{lccccccccccccl}\nonumber{}\mbox{maximize} & x_{0} &  &  & + & x_{2} & + &  & - & x_{4} & + & x_{5} &  & \\\nonumber{}\mbox{subject to} & -x_{0} &  &  &  &  & + & x_{3} &  &  &  &  & = & 1,\\\nonumber{} &  &  &  &  & x_{2} & - & x_{3} & + & x_{4} & + & x_{5} & = & -2,\\\nonumber{} & x_{0} & - & x_{1} &  &  &  &  & - & x_{4} & - & x_{5} & = & 0,\\\nonumber{} &  &  & x_{1} & - & x_{2} & + &  &  &  &  &  & = & 0,\end{array}\end{math} (6.1.2)

having the bounds [[MathCmd 132]].

The corresponding graph [[MathCmd 119]] is displayed in fig.6.1.

Figure 6.1: Simple network.

6.1.1.1. Source code

In this section we will show how to solve (6.1.2) with the network optimizer.

The Python program included below, which solves this problem, is distributed with MOSEK and can be found in the directory

  mosek\6\tools\examples\python

# # Copyright: Copyright (c) 1998-2011 MOSEK ApS, Denmark. All rights reserved. # # File: network1.py # # Demonstrates a simple use of the network optimizer. # # Purpose: 1. Specify data for a network. # 2. Solve the network problem with the network optimizer. # import sys import mosek from numpy import array, float, zeros # Define a stream printer to grab output from MOSEK def streamprinter(text): sys.stdout.write(text) sys.stdout.flush() # Since the actual value of Infinity is ignored, we define it solely # for symbolic purposes: inf = 0.0 # Make mosek environment. env = mosek.Env () # Initialize the environment. env.init () # Attach a printer to the environment env.set_Stream (mosek.streamtype.log, streamprinter) numcon = 4 numvar = 6 # Specify network graph data. netfrom = array ([0,2,3,1,1,1]) netto = array ([2,3,1,0,2,2]) # Specify arc cost. cc = zeros (4, float) cx = array ([1.0,0.0,1.0,0.0,-1.0,1.0]) # Specify boundkeys. bkc = [mosek.boundkey.fx]*4 bkx = [mosek.boundkey.lo]*6 # Specify bounds. blc = array ([1.0,1.0,-2.0,0.0]) buc = array ([1.0,1.0,-2.0,0.0]) blx = zeros (6, float) bux = array ([inf,inf,inf,inf,inf,inf]) # Specify zero primal solution. xc = zeros (4, float) xx = zeros (6, float) # Specify zero dual solution. y = zeros (4, float) slc = zeros (4, float) suc = zeros (4, float) slx = zeros (6, float) sux = zeros (6, float) # Specify status keys. skc = [mosek.stakey.unk]*4 skx = [mosek.stakey.unk]*6 # Create a task object linked with the environment env. dummytask = env.Task (numcon,numvar) # Set the problem to be maximized dummytask.putobjsense (mosek.objsense.maximize) # Solve the network problem prosta,solsta = dummytask.netoptimize( cc, cx, bkc, blc, buc, bkx, blx, bux, netfrom, netto, 0, skc, skx, xc, xx, y, slc, suc, slx, sux) if solsta == mosek.solsta.optimal : print "Network problem is optimal" print "Primal solution is :" for i in range(0,numcon) : print "xc[%d] = %-16.10e" % (i,xc[i]) for j in range(0,numvar) : print "Arc(%d,%d) -> xx[%d] = %-16.10e" % (netfrom[j],netto[j],j,xx[j]) elif solsta == mosek.solsta.prim_infeas_cer : print "Network problem is primal infeasible" elif solsta == mosek.solsta.dual_infeas_cer : print "Network problem is dual infeasible" else : print "Network problem solsta : %s" % solsta

6.1.1.2. Example code comments

There are a few important differences between the linear network optimization example in section 6.1.1.1 and the general linear optimization problem in section 5.2.

  • MOSEK allows that network problems can be inputted and optimized using one function call to the function Task.netoptimize. This is more efficient and uses less memory than a call to the standard optimizer.
  • Since we know that each column of matrix A has two non-zeroes, it can be stored in two arrays, from and to, specifying the origin and destination of the arcs (variables), see graph in fig.fig-network.
  • The solution is written directly to skc, skx, xc, xx, y, slc, suc, slx and sux by Task.netoptimize.

6.2. Embedded network flow problems

Often problems contains both large parts with network structure and some non-network constraints or variables — such problems are said to have embedded network structure.

A linear optimization with embedded network structure problem can be written as :

\begin{math}\nonumber{}\begin{array}{lccccl}\nonumber{}\mbox{minimize} &  &  & c^{T}x+c^{f} &  & \\\nonumber{}\mbox{subject to} & l^{c}_{N} & \leq{} & Nx & \leq{} & u^{c}_{N},\\\nonumber{} & l^{c} & \leq{} & Ax & \leq{} & u^{c},\\\nonumber{} & l^{x} & \leq{} & x & \leq{} & u^{x},\end{array}\end{math} (6.2.1)

Where the constraints

\begin{math}\nonumber{}\begin{array}{lcccl}\nonumber{}l^{c}_{N} & \leq{} & Nx & \leq{} & u^{c}_{N}\end{array}\end{math} (6.2.2)

defines a network as explained in section 6.1, and the constraints

\begin{math}\nonumber{}\begin{array}{lcccl}\nonumber{}l^{c} & \leq{} & Ax & \leq{} & u^{c}\end{array}\end{math} (6.2.3)

defines the general non-network linear constraints. As an example consider the small linear optimization problem

\begin{math}\nonumber{}\begin{array}{lccccccccccccl}\nonumber{}\mbox{maximize} & -x_{0} &  &  & + & x_{2} &  &  & - & x_{4} & + & x_{5} &  & \\\nonumber{}\mbox{subject to} & 0.50x_{0} &  &  &  &  & + & 0.50x_{3} &  &  &  &  & = & 0.5,\\\nonumber{} &  &  &  &  & 0.50x_{2} & - & 0.50x_{3} & + & 0.50x_{4} & + & 0.50x_{5} & = & -1,\\\nonumber{} & -0.25x_{0} & + & -2.50x_{1} & + &  &  &  & - & 0.25x_{4} & - & 0.25x_{5} & = & 0,\\\nonumber{} &  &  & 2.50x_{1} & - & 0.25x_{2} &  &  &  &  &  &  & = & 0,\\\nonumber{} &  & - & x_{1} & + & x_{2} & + & x_{3} &  &  & + & x_{5} & \geq{} & 6,\end{array}\end{math} (6.2.4)

with the bounds

\begin{displaymath}\nonumber{}-\infty \leq{}x_{0}\leq{}0,0\leq{}x_{j}\leq{}\infty \mbox{ for }j=1\ldots 5.\end{displaymath}

Recalling the network flow problem structural properties from section 6.1.1, each variable should appear in exactly two constraints with coefficients of either [[MathCmd 129]] or [[MathCmd 130]].

At first glance it does not seem to contain any network structure, but if we scale constraints 1-4 by respectively 2.0, 2.0, 4.0, 4.0 and columns 1-2 by -1.0, 0.1 we get the following problem :

\begin{math}\nonumber{}\begin{array}{lccccccccccccl}\nonumber{}\mbox{maximize} & x_{0} &  &  & + & x_{2} & + &  & - & x_{4} & + & x_{5} &  & \\\nonumber{}\mbox{subject to} & -x_{0} &  &  &  &  & + & x_{3} &  &  &  &  & = & 1,\\\nonumber{} &  &  &  &  & x_{2} & - & x_{3} & + & x_{4} & + & x_{5} & = & -2,\\\nonumber{} & x_{0} & - & x_{1} &  &  &  &  & - & x_{4} & - & x_{5} & = & 0,\\\nonumber{} &  &  & x_{1} & - & x_{2} & + &  &  &  &  &  & = & 0,\\\nonumber{} &  &  & x_{1} & + & x_{2} & + & x_{3} &  &  & + & x_{5} & \geq{} & 6,\end{array}\end{math} (6.2.5)

with the bounds

\begin{displaymath}\nonumber{}0\leq{}x_{j}\leq{}\infty \mbox{ for }j=0\ldots 5.\end{displaymath}

This corresponds to the network flow problem in section 6.1.1 plus one extra non-network constraint. We cannot use the network optimizer directly on the above problem since the last constraint destroys the network property. Finding the largest possible network structure in a linear optimization problem is computationally difficult, so MOSEK offers a heuristic Task.netextraction that attempts to find suitable scaling factors maximizing numbers of network constraints and variables. Assuming that the embedded network structure is dominant and the problem has few non-network constraints, we can exploit this structure and potentially speed up the optimization. Since the network constraints can be handled efficiently by the specialized network optimizer, the following idea is used:

An embedded network can be exploited by this scheme in two ways:

The first method is more difficult than the second, but also offers much more flexibility. In 6.2.1 the first method is demonstrated by a code example below. For further details on exploiting embedded network structure in the standard simplex optimizer, see section 8.3.1.

6.2.1. Example: Exploit embedded network flow structure in the simplex optimizer

MOSEK is distributed with some network examples which can be found in the directory

  mosek\6\tools\examples 

The example given in this section demonstrates how to extract and optimize embedded network structure in a arbitrary linear optimization problem. The following idea is used

  • Read an arbitrary linear optimization problem into a task.
  • Use the Task.netextraction function to extract embedded network structure.
  • Optimize the network problem using the Task.netoptimize function.

# # Copyright: Copyright (c) 1998-2011 MOSEK ApS, Denmark. All rights reserved. # # File: network2.py # # Demonstrates a simple use of network structure in a model. # # Purpose: 1. Read an optimization problem from an # user specified MPS file. # 2. Extract the embedded network (if any ). # 3. Solve the embedded network with the network optimizer. # # Note that the general simplex optimizer called though MSK_optimize can also extract # embedded network and solve it with the network optimizer. The direct call to the # network optimizer, which is demonstrated here, is offered as an option to save # memory and overhead for solving either many or large network problems. # import sys import mosek from numpy import resize,array, float, zeros #from mosek.array import * # Define a stream printer to grab output from MOSEK def streamprinter(text): sys.stdout.write(text) sys.stdout.flush() # Since the actual value of Infinity is ignores, we define it solely # for symbolic purposes: inf = 0.0 if len(sys.argv) != 2: print "Wrong arguments. The syntax is:" print " network2 inputfile" else: # Make mosek environment. env = mosek.Env () # Initialize the environment. env.init () # Attach a printer to the environment env.set_Stream (mosek.streamtype.log, streamprinter) # since we don't know the size of the problem. numcon = [] numvar = [] netnumcon = [] netnumvar = [] task = env.Task (0,0) task.readdata (sys.argv[1]) numcon = task.getnumcon () numvar = task.getnumvar () # Specify network graph data. netfrom = zeros (numvar, int) netto = zeros (numvar, int) # Specify arc cost. cc = zeros (numcon, float) cx = zeros (numvar, float) # Specify boundkeys. bkc = [ mosek.boundkey.fx ] * numcon bkx = [ mosek.boundkey.fx ] * numvar) # Specify bounds. blc = zeros (numcon, float) buc = zeros (numcon, float) blx = zeros (numvar, float) bux = zeros (numvar, float) # Specify data for extracted network. scalcon = zeros (numcon, float) scalvar = zeros (numvar, float) netcon = zeros (numcon, int) netvar = zeros (numvar, int) # Extract embedded network netnumcon,netnumvar = task.netextraction( netcon, netvar, scalcon, scalvar, cx, bkc, blc, buc, bkx, blx, bux, netfrom, netto) # Create a dummy task object linked with the environment env. dummytask = env.Task (netnumcon,netnumvar) # Array length for netoptimize must match netnumcon and netnumvar # Resize network graph data. netfrom = resize (netfrom,netnumvar) netto = resize (netto,netnumvar) # Resize arc cost. cc = resize (cc,netnumcon) cx = resize (cx,netnumvar) # Resize boundkeys. bkc = [ mosek.boundkey.fx ] * netnumcon bkx = [ mosek.boundkey.fx ] * netnumvar # Resize bounds. blc = resize (blc,netnumcon) buc = resize (buc,netnumcon) blx = resize (blx,netnumvar) bux = resize (bux,netnumvar) # Specify zero primal solution. xc = zeros (netnumcon, float) xx = zeros (netnumvar, float) # Specify zero dual solution. y = zeros (netnumcon, float) slc = zeros (netnumcon, float) suc = zeros (netnumcon, float) slx = zeros (netnumvar, float) sux = zeros (netnumvar, float) # Specify status keys. skc = [ mosek.stakey.unk ] * netnumcon skx = [ mosek.stakey.unk ] * netnumvar # Specify problem and solution status. prosta = [] solsta = [] # Solve the network problem prosta,solsta = dummytask.netoptimize( cc, cx, bkc, blc, buc, bkx, blx, bux, netfrom, netto, 0, skc, skx, xc, xx, y, slc, suc, slx, sux) print "Original problem size : numcon : %d numvar : %d" % (numcon,numvar) print "Embedded network size : numcon : %d numvar : %d" % (netnumcon,netnumvar) if solsta == mosek.solsta.optimal : print "Network problem is optimal" elif solsta == mosek.solsta.prim_infeas_cer : print "Network problem is primal infeasible" elif solsta == mosek.solsta.dual_infeas_cer : print "Network problem is dual infeasible" else : print "Network problem solsta : %s" % solsta

In the above example we only optimize the embedded network problem. We still need to use the found network solution as a hot-start for the simplex optimizer and solve the original problem. This involves unscaling the network solution back to same unit measure as the original problem. In the example

  mosek\6\tools\examples\python\network3.py

we show how to convert the network solution into a valid hot-start for the simplex optimizer.

Wed Feb 29 16:13:01 2012