13. Case Studies


13.1. The traveling salesman problem

The Travelling Salesman Problem (TSP) is the problem of finding the shortest cyclic tour between a set of cities, visiting each city exactly once. This can be formulated using mixed integer programming.

When solving mixed integer optimization problems it is important to use a strong formulation of the problem, otherwise MOSEK may spend a very long time solving the optimization problem. This is not only true for MOSEK but for the branch-and-bound based solution method too.

The approach explored in this section is an implementation of the approach discussed in the article “Teaching integer programming formulations using the Traveling Salesman Problem” by Gábor Pataki [13].

13.1.1. The TSP formulations

Given a set of nodes we want to find the shortest tour (a directed cycle containing all nodes) in a complete directed graph. We use the variables [[MathCmd 608]] to indicate whether the arc (i,j) is included in the tour.

The core of the formulation is

\begin{math}\nonumber{}\begin{array}{lll}\nonumber{}\mbox{minimize} & \sum \limits _{{i,j}}c_{{ij}}x_{{ij}} & \\\nonumber{}\mbox{subject to} & \sum \limits _{i}x_{{ij}}=1 & \forall j,\\\nonumber{} & \sum \limits _{j}x_{{ij}}=1 & \forall i,\\\nonumber{} & 0\leq{}x_{{ij}}\leq{}1, & x_{{ij}}\,\mbox{integer}.\end{array}\end{math} (13.1.1)

These constraints are called the assignment constraints. The assignment constraints, however, do not constitute the entire formulation as groups of disjoint cycles, called subtours, as well as the complete tours are feasible.

To exclude the subtours two sets of constraints are considered.

The MTZ formulation

The MTZ (Miller-Tucker-Zemlin) formulation of the TSP includes the following constraints

\begin{math}\nonumber{}\begin{array}{cl}\nonumber{}u_{1}=1, & \\\nonumber{}2\leq{}u_{i}\leq{}n & \forall i\neq{}1,\\\nonumber{}u_{i}-u_{j}+1\leq{}(n-1)(1-x_{{ij}}) & \forall i\neq{}1,j\neq{}1.\end{array}\end{math} (13.1.2)

The idea of this formulation is to assign the numbers 1 through n to the nodes with the extra variables [[MathCmd 802]] so that this numbering corresponds to the order of the nodes in the tour. It is obvious that this excludes subtours, as a subtour excluding the node 1 cannot have a feasible assignment of the corresponding [[MathCmd 802]] variables.

The subtour formulation

An alternative approach is simply to take any potential subtour, i.e. any true subset of nodes, and declare that it is illegal.

\begin{math}\nonumber{}\sum _{{i,j\in{}S}}x_{{ij}}\leq{}|S|-1\quad{}(S\not\subseteq{}V,|S|>1)\end{math} (13.1.3)

As the subtour inequality for [[MathCmd 805]] is a linear combination of the inequality for S and the assignment constraints, it is sufficient to use the subtour inequalities with S having size n/2 at most. Note that this formulation has the disadvantage of being exponential in size.

The MTZ formulation of the TSP is a very weak formulation so we will try to strengthen it by adding some of the subtour constraints from the stronger subtour formulation and then compare the solution times. For each problem we will try to identify some of the most relevant subtour constraints by solving the relaxed IP without the MTZ constraints and then choosing some of the violated subtour inequalities, corresponding to the subtours in the solution. The complete algorithm in pseudo-code is (the complete C implementation is included below in Section 13.1.3):

  1. Let the IP formulation consist of the assignment constraints (13.1.1) only.
  2. for k=1 to [[MathCmd 806]]

    • Solve the IP over the current formulation. Assume that the optimal solution consists of r subtours [[MathCmd 807]].
    • If r=1, stop; the solution is optimal to the TSP. Otherwise, add to the formulation 1000 subtour constraints at most, in which S is the union of several [[MathCmd 808]] sets and [[MathCmd 809]].
  3. Add the MTZ arc constraints to the formulation, and solve the IP to optimality.

Each round we add 1000 constraints at most as the number of violated subtour inequalities is exponential in r.

Setting [[MathCmd 806]] equal to 0, 1, and 2, we obtain three formulations of increasing strength which we solve in 3.

13.1.2. Comparing formulations

We have tested this method on six TSP instances from the TSPLIB library which can be found at

http://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/

The time and number of B&B nodes for each of the three formulations are recorded in Table 13.1. The entry “***” means that the problem was unsolvable within a time window of 5000 seconds. The time spent solving the relaxed IPs and identifying subtour constraints was negligible.

Number of rounds Zero rounds One round Two rounds
Problem name Time B&B nodes Time B&B nodes Time B&B nodes
bays29 658 53809 85 1715 39 2739
berlin52 553 7944 56 198 10 2
br17 *** *** 1 13 1 1
ft70 *** *** 17 3 16 5
ftv33 23 1882 8 2 9 1
ftv55 864 12494 138 4853 53 2515
Table 13.1: Solving TSP using increasingly stronger formulations.

Not surprisingly a stronger formulation means shorter solution time (with a few exceptions where the second round of strengthening seemingly is superfluous), but it is worth noting the magnitude of the decrease in solution time arising from stronger formulations.

Therefore, it is often worthwhile to consider whether one can strengthen a given formulation when solving a mixed integer optimization problem.

13.1.3. Example code

The following example is included in the distribution in the file msktsp.c.

/* Copyright: Copyright (c) 1998-2011 MOSEK ApS, Denmark. All rights reserved. File: msktsp.c Purpose: Demonstrates the difference between weak and strong formulations when solving MIP's. */ #include <stdio.h> #include <string.h> #include <math.h> #include <assert.h> #include "mosek.h" #define MAXCUTROUNDS 2 #define MAXADDPERROUND 1000 static void MSKAPI printstr(void *handle, char str[]) { printf("MOSEK: %s",str); } /* printstr */ /* conversion from n x n tsp city matrix indices to array index */ #define IJ(i,j) (n*(i)+(j)) /* mallocs and returns costmatrix, returns number of cities in ncities */ int* readtspfromfile(char* filename, int* ncities) { FILE *tspfile; char sbuf[21]; tspfile = fopen(filename,"r"); if (!tspfile) return NULL; do { if (1 != fscanf(tspfile,"%20s ",sbuf)) return NULL; } while (strncmp(sbuf,"DIMENSION",9) != 0); if (1 != fscanf(tspfile,"%d ",ncities)) return NULL; do { if (1 != fscanf(tspfile,"%20s ",sbuf)) return NULL; } while (strncmp(sbuf,"EDGE_WEIGHT_TYPE",16) != 0); if (1 != fscanf(tspfile,"%20s ",sbuf)) return NULL; if (strcmp(sbuf,"EXPLICIT") == 0) { do { if (1 != fscanf(tspfile,"%20s ",sbuf)) return NULL; } while (strncmp(sbuf,"EDGE_WEIGHT_FORMAT",18) != 0); if (1 != fscanf(tspfile,"%20s ",sbuf)) return NULL; if (strcmp(sbuf,"FULL_MATRIX") == 0) { int* cost; int ij, n2; do { if (1 != fscanf(tspfile,"%20s ",sbuf)) return NULL; } while (strncmp(sbuf,"EDGE_WEIGHT_SECTION",19) != 0); n2 = *ncities; n2 *= n2; cost = (int*) malloc(n2*sizeof(int)); assert(cost); for (ij = 0; ij<n2; ij++) { if (1 != fscanf(tspfile,"%d ",&cost[ij])) { free(cost); return NULL; } } return cost; } else if (strcmp(sbuf,"LOWER_DIAG_ROW") == 0) { int* cost; int i, j, n; do { if (1 != fscanf(tspfile,"%20s ",sbuf)) return NULL; } while (strncmp(sbuf,"EDGE_WEIGHT_SECTION",19) != 0); n = *ncities; cost = (int*) malloc(n*n*sizeof(int)); assert(cost); for (i=0; i<n; i++) for (j=0; j<=i; j++) { int c; if (1 != fscanf(tspfile,"%d ",&c)) { free(cost); return NULL; } cost[IJ(i,j)] = c; cost[IJ(j,i)] = c; } return cost; } else { printf("Format not supported\n"); return NULL; } } else if (strcmp(sbuf,"EUC_2D") == 0) { int* cost; double *xcoord, *ycoord; int i, j, n; do { if (1 != fscanf(tspfile,"%20s ",sbuf)) return NULL; } while (strncmp(sbuf,"NODE_COORD_SECTION",18) != 0); n = *ncities; xcoord = (double*) malloc(n*sizeof(double)); ycoord = (double*) malloc(n*sizeof(double)); cost = (int*) malloc(n*n*sizeof(int)); assert(xcoord); assert(ycoord); assert(cost); for (i = 0; i<n; i++) { int dummy; if (3 != fscanf(tspfile,"%d %lf %lf ",&dummy,&xcoord[i],&ycoord[i])) { free(cost); return NULL; } } for (i = 0; i<n; i++) for (j=0; j<n; j++) { double xd = xcoord[i] - xcoord[j]; double yd = ycoord[i] - ycoord[j]; cost[IJ(i,j)] = (int) (0.5 + sqrt(xd*xd + yd*yd)); } return cost; } else { printf("E_W_Type not supported\n"); return NULL; } } /* readtspfromfile */ /* add the x_ij variables */ void add_vars(MSKtask_t task, int n) { MSKrescodee r; int ij; int n2 = n*n; r = MSK_append(task,MSK_ACC_VAR,n2); assert(r==MSK_RES_OK); for(ij=0; ij<n2; ++ij) { r = MSK_putbound(task,MSK_ACC_VAR,ij,MSK_BK_RA,0,1); assert(r==MSK_RES_OK); r = MSK_putvartype(task,ij,MSK_VAR_TYPE_INT); assert(r==MSK_RES_OK); } for (ij=0; ij<n; ij++) { r = MSK_putbound(task,MSK_ACC_VAR,IJ(ij,ij),MSK_BK_FX,0,0); assert(r==MSK_RES_OK); } } /* add_vars */ /* adds the tsp objective function and frees cost */ void add_objective_function(MSKtask_t task, int n, int* cost) { MSKrescodee r; int ij; int n2 = n*n; r = MSK_putcfix(task,0.0); assert(r==MSK_RES_OK); for(ij=0; ij<n2; ++ij) { r = MSK_putcj(task,ij,cost[ij]); assert(r==MSK_RES_OK); } free(cost); } /* add_objective_function */ /* adds the tsp assignment constraints */ void add_assignment_constraints(MSKtask_t task, int n) { MSKrescodee r; int i, j; double* aval; int *asub; aval = (double*) malloc(n*sizeof(double)); assert(aval); asub = (int*) malloc(n*sizeof(int)); assert(asub); for (i=0; i<n; i++) aval[i] = 1; r = MSK_append(task,MSK_ACC_CON,n*2); assert(r==MSK_RES_OK); /* Constraint 0--(n-1) is \sum_j x_{ij} = 1 */ for (i=0; i<n; i++) { r = MSK_putbound(task,MSK_ACC_CON,i,MSK_BK_FX,1,1); assert(r==MSK_RES_OK); for (j=0; j<n; j++) asub[j] = IJ(i,j); r = MSK_putavec(task,MSK_ACC_CON,i,n,asub,aval); assert(r==MSK_RES_OK); } /* Constraint n--(2n-1) is \sum_i x_{ij} = 1 */ for (j=0; j<n; j++) { r = MSK_putbound(task,MSK_ACC_CON,j+n,MSK_BK_FX,1,1); assert(r==MSK_RES_OK); for (i=0; i<n; i++) asub[i] = IJ(i,j); r = MSK_putavec(task,MSK_ACC_CON,j+n,n,asub,aval); assert(r==MSK_RES_OK); } free(aval); free(asub); } /* add_assignment_constraints */ /* adds the Miller-Tucker-Zemlin arc constraints */ void add_MTZ_arc_constraints(MSKtask_t task, int n) { MSKrescodee r; int varidx, conidx, i, j; r = MSK_getnumvar(task,&varidx); assert(r==MSK_RES_OK); r = MSK_getnumcon(task,&conidx); assert(r==MSK_RES_OK); /* add the vars u_k for k=1..(n-1) getting index * from varidx to varidx+n-2 */ r = MSK_append(task,MSK_ACC_VAR,n-1); assert(r==MSK_RES_OK); for(i=varidx; i<varidx+n-1; ++i) { /* set bound: 2 <= u_k <= n, k=1..(n-1) */ r = MSK_putbound(task,MSK_ACC_VAR,i,MSK_BK_RA,2,n); assert(r==MSK_RES_OK); } /* add the (n-1)^2 constraints: * u_i - u_j + 1 <= (n - 1)(1 - x_ij) or equivalently * u_i - u_j + (n - 1)x_ij <= n - 2, for i,j != 0 */ r = MSK_append(task,MSK_ACC_CON,(n-1)*(n-1)); assert(r==MSK_RES_OK); for (i=1; i<n; i++) for (j=1; j<n; j++) { double aval[3]; int asub[3]; aval[0] = 1; aval[1] = -1; aval[2] = n-1; asub[0] = varidx + i - 1; /* u_i */ asub[1] = varidx + j - 1; /* u_j */ asub[2] = IJ(i,j); /* x_ij */ r = MSK_putbound(task,MSK_ACC_CON,conidx,MSK_BK_UP,-MSK_INFINITY,n-2); assert(r==MSK_RES_OK); r = MSK_putavec(task,MSK_ACC_CON,conidx,3,asub,aval); assert(r==MSK_RES_OK); conidx++; } } /* add_MTZ_arc_constraints */ /* construct the list of cities in the chosen subtours */ int* subtourstolist(MSKtask_t task, int n, int nextnode[], int subtour[], int chosen[], int k, int* size) { int ncities, i, j; int *cities; cities = (int*) malloc(n*sizeof(int)); assert(cities); ncities = 0; for (i=0; i<k; i++) { int subtourstart = subtour[chosen[i]]; j = subtourstart; do { cities[ncities] = j; ncities++; j = nextnode[j]; } while (j != subtourstart); } *size = ncities; return cities; } /* subtourstolist */ /* adds the subtour constraint given by the list cities S: * \sum_{i,j \in S} x_{ij} \leq |S|-1 */ void addcut(MSKtask_t task, int n, int citylist[], int size) { MSKrescodee r; int i, j, asubidx, conidx; double* aval; int *asub; int size2 = size*size; aval = (double*) malloc(size2*sizeof(double)); assert(aval); asub = (int*) malloc(size2*sizeof(int)); assert(asub); for (i=0; i<size2; i++) aval[i] = 1; r = MSK_getnumcon(task,&conidx); assert(r==MSK_RES_OK); r = MSK_append(task,MSK_ACC_CON,1); assert(r==MSK_RES_OK); r = MSK_putbound(task,MSK_ACC_CON,conidx,MSK_BK_UP,-MSK_INFINITY,size-1); assert(r==MSK_RES_OK); asubidx = 0; for (i=0; i<size; i++) for (j=0; j<size; j++) { asub[asubidx] = IJ(citylist[i],citylist[j]); asubidx++; } r = MSK_putavec(task,MSK_ACC_CON,conidx,size2,asub,aval); assert(r==MSK_RES_OK); free(aval); free(asub); } /* addcut */ /* identifies subtours and adds a number of violated cuts */ void addcuts(MSKtask_t task, int n, int maxcuts, int* nsubtours, int* ncuts) { MSKrescodee r; int i, j, k; int n2 = n*n; double *xx; int *nextnode, *visited, *subtour, *chosen; int nsubt = 0; xx = (double*) malloc(n2*sizeof(double)); nextnode = (int*) malloc(n*sizeof(int)); assert(xx); assert(nextnode); r = MSK_getsolutionslice(task,MSK_SOL_ITG,MSK_SOL_ITEM_XX,0,n2,xx); assert(r==MSK_RES_OK); /* convert matrix representation of graph (xx) to * adjacency(-list) (nextnode) */ for (i=0; i<n; i++) for (j=0; j<n; j++) { if (xx[IJ(i,j)]>0.5) /* i.e. x_ij = 1 */ nextnode[i] = j; } free(xx); xx = NULL; visited = (int*) calloc(n,sizeof(int)); /* visited is initialized to 0 */ subtour = (int*) malloc(n*sizeof(int)); assert(visited); assert(subtour); /* identify subtours; keep count in nsubt, save starting * pointers in subtour[0..(nsubt-1)] */ for (i=0; i<n; i++) if (!visited[i]) /* find an unvisited node; * this starts a new subtour */ { subtour[nsubt] = i; nsubt++; j = i; do { assert(!visited[j]); visited[j] = 1; j = nextnode[j]; } while (j!=i); } free(visited); visited = NULL; *nsubtours = nsubt; *ncuts = 0; chosen = (int*) malloc(nsubt*sizeof(int)); /* list of chosen subtours */ for (k=1; k<=nsubt; k++) /* choose k of nsubt subtours */ { int nchosen = 1; chosen[0] = nsubt - 1; while (*ncuts < maxcuts) { if (nchosen == k) { int *citylist; int size; citylist = subtourstolist(task,n,nextnode,subtour, chosen,k,&size); if (size <= n/2) /* add only subtour constraints * of size n/2 or less */ { addcut(task,n,citylist,size); (*ncuts)++; } free(citylist); j=0; while (j<k && chosen[k - 1 - j] == j) j++; if (k==j) break; /* all k-size subsets done */ nchosen = k - j; chosen[nchosen - 1]--; } else /* 0 < nchosen < k */ { chosen[nchosen] = chosen[nchosen - 1] - 1; nchosen++; } } } free(nextnode); free(subtour); free(chosen); } /* addcuts */ int main(int argc, char *argv[]) { int *cost; /* tsp cost matrix */ int n; /* number of cities */ MSKenv_t env; /* Mosek environment */ MSKtask_t task; /* Mosek task */ MSKrescodee r; /* Mosek return code */ double ObjVal; /* Value of the objective function */ int maxrounds; /* number of cutting rounds */ int maxcuts; /* maximum number of cuts added per round */ int k; int nsubtours, ncuts; double t; double cuttime = 0; if (argc < 2) { printf("Usage: ./tsp filename.tsp [rounds] [maxcuts]\n\n" "rounds is the maximum number of cutting rounds (default = %d)\n" "maxcuts is the maximum number of cuts added per round " "(default = %d)\n", MAXCUTROUNDS,MAXADDPERROUND); return 1; } maxrounds = MAXCUTROUNDS; if (argc >= 3) maxrounds = atoi(argv[2]); maxcuts = MAXADDPERROUND; if (argc >= 4) maxcuts = atoi(argv[3]); cost = readtspfromfile(argv[1],&n); if (!cost) { printf("Bad tsp file\n"); return 1; } r = MSK_makeenv(&env,NULL,NULL,NULL,NULL); assert(r==MSK_RES_OK); MSK_linkfunctoenvstream(env,MSK_STREAM_LOG,NULL,printstr); r = MSK_initenv(env); assert(r==MSK_RES_OK); r = MSK_makeemptytask(env,&task); assert(r==MSK_RES_OK); MSK_linkfunctotaskstream(task,MSK_STREAM_LOG,NULL,printstr); add_vars(task,n); add_objective_function(task,n,cost); add_assignment_constraints(task,n); nsubtours = 2; for (k=0; k<maxrounds; k++) { r = MSK_optimize(task); assert(r==MSK_RES_OK); r = MSK_getprimalobj(task,MSK_SOL_ITG,&ObjVal); assert(r==MSK_RES_OK); MSK_getdouinf(task,MSK_DINF_OPTIMIZER_TIME,&t); cuttime += t; addcuts(task,n,maxcuts,&nsubtours,&ncuts); printf("\n" "Round: %d\n" "ObjValue: %e\n" "Number of subtours: %d\n" "Number of cuts added: %d\n\n",k+1,ObjVal,nsubtours,ncuts); if (nsubtours == 1) break; /* problem solved! */ } t = 0; if (nsubtours > 1) { printf("Adding MTZ arc constraints\n\n"); add_MTZ_arc_constraints(task,n); r = MSK_optimize(task); assert(r==MSK_RES_OK); r = MSK_getprimalobj(task,MSK_SOL_ITG,&ObjVal); assert(r==MSK_RES_OK); MSK_getdouinf(task,MSK_DINF_OPTIMIZER_TIME,&t); } printf("\n" "Done solving.\n" "Time spent cutting: %.2f\n" "Total time spent: %.2f\n" "ObjValue: %e\n",cuttime,cuttime+t,ObjVal); MSK_deletetask(&task); MSK_deleteenv(&env); return 0; } /* main */

13.2. Geometric (posynomial) optimization

13.2.1. The problem

A geometric optimization problem can be stated as follows

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{k\in{}J_{0}}}c_{k}\prod \limits _{{j=0}}^{{n-1}}t_{j}^{{a_{{kj}}}} &  &  & \\\nonumber{}\mbox{subject to} & \sum \limits _{{k\in{}J_{i}}}c_{k}\prod \limits _{{j=0}}^{{n-1}}t_{j}^{{a_{{kj}}}} & \leq{} & 1, & i=1,\ldots ,m,\\\nonumber{} & t>0, &  &  &\end{array}\end{math} (13.2.1)

where it is assumed that

\begin{displaymath}\nonumber{}\cup _{{k=0}}^{m}J_{k}=\lbrace{}1,\ldots ,T\rbrace{}\end{displaymath}

and if [[MathCmd 169]], then

\begin{displaymath}\nonumber{}J_{i}\cap J_{j}=\emptyset .\end{displaymath}

Hence, A is a [[MathCmd 815]] matrix and c is a vector of length T. Given [[MathCmd 816]] then

\begin{displaymath}\nonumber{}c_{k}\prod \limits _{{j=0}}^{{n-1}}t_{j}^{{a_{{kj}}}}\end{displaymath}

is called a monomial . A sum of monomials i.e.

\begin{displaymath}\nonumber{}\sum \limits _{{k\in{}J_{i}}}c_{k}\prod \limits _{{j=0}}^{{n-1}}t_{j}^{{a_{{kj}}}}\end{displaymath}

is called a posynomial . In general, the problem (13.2.1) is very hard to solve. However, the posynomial case where it is required that

\begin{displaymath}\nonumber{}c>0\end{displaymath}

is relatively easy. The reason is that using a simple variable transformation a convex optimization problem can be obtained. Indeed using the variable transformation

\begin{math}\nonumber{}t_{j}=e^{{x_{j}}}\end{math} (13.2.2)

we obtain the problem

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{k\in{}J_{0}}}c_{k}e^{{\sum \limits _{{j=0}}^{{n-1}}a_{{kj}}x_{j}}} &  &  & \\\nonumber{}\mbox{subject to} & \sum \limits _{{k\in{}J_{i}}}c_{k}e^{{\sum \limits _{{j=0}}^{{n-1}}a_{{kj}}x_{j}}} & \leq{} & 1, & i=1,\ldots ,m,\end{array}\end{math} (13.2.3)

which is a convex optimization problem that can be solved using MOSEK. We will call

\begin{displaymath}\nonumber{}c_{t}e^{{\left(\sum \limits _{{j=0}}^{{n-1}}a_{{tj}}x_{j}\right)}}=e^{{\left(\log (c_{t})+\sum \limits _{{j=0}}^{{n-1}}a_{{tj}}x_{j}\right)}}\end{displaymath}

a term and hence the number of terms is T.

As stated, the problem (13.2.3) is non-separable. However, using

\begin{displaymath}\nonumber{}v_{t}=\log (c_{t})+\sum \limits _{{j=0}}^{{n-1}}a_{{tj}}x_{j}\end{displaymath}

we obtain the separable problem

\begin{math}\nonumber{}\begin{array}{lccll}\nonumber{}\mbox{minimize} & \sum \limits _{{t\in{}J_{0}}}e^{{v_{t}}} &  &  & \\\nonumber{}\mbox{subject to} & \sum \limits _{{t\in{}J_{i}}}e^{{v_{t}}} & \leq{} & 1, & i=1,\ldots ,m,\\\nonumber{} & \sum \limits _{{j=0}}^{{n-1}}a_{{tj}}x_{j}-v_{t} & = & -\log (c_{t}), & t=0,\ldots ,T,\end{array}\end{math} (13.2.4)

which is a separable convex optimization problem.

A warning about this approach is that the exponential function [[MathCmd 153]] is only numerically well-defined for values of x in a small interval around 0 since [[MathCmd 153]] grows very rapidly as x becomes larger. Therefore numerical problems may arise when solving the problem on this form.

13.2.2. Applications

A large number of practical applications, particularly in electrical circuit design, can be cast as a geometric optimization problem. We will not review these applications here but rather refer the reader to [14] and the references therein.

13.2.3. Modeling tricks

A lot of tricks that can be used for modeling posynomial optimization problems are described in [14]. Therefore, in this section we cover only one important case.

13.2.3.1. Equalities

In general, equalities are not allowed in (13.2.1), i.e.

\begin{displaymath}\nonumber{}\sum \limits _{{k\in{}J_{i}}}c_{k}\prod \limits _{{j=0}}^{{n-1}}t_{j}^{{a_{{kj}}}}=1\end{displaymath}

is not allowed. However, a monomial equality is not a problem. Indeed consider the example

\begin{displaymath}\nonumber{}xyz^{{-1}}=1\end{displaymath}

of a monomial equality. The equality is identical to

\begin{displaymath}\nonumber{}1\leq{}xyz^{{-1}}\leq{}1\end{displaymath}

which in turn is identical to the two inequalities

\begin{displaymath}\nonumber{}\begin{array}{lclcl}\nonumber{}xyz^{{-1}} &  &  & \leq{} & 1,\\\nonumber{}\frac{1}{xyz^{{-1}}} & = & x^{{-1}}y^{{-1}}z & \leq{} & 1.\end{array}\end{displaymath}

Hence, it is possible to model a monomial equality using two inequalities.

13.2.4. Problematic formulations

Certain formulations of geometric optimization problems may cause problems for the algorithms implemented in MOSEK. Basically there are two kinds of problems that may occur:

  • The solution vector is finite, but an optimal objective value can only be a approximated.
  • The optimal objective value is finite but implies that a variable in the solution is infinite.

13.2.4.1. Finite unattainable solution

The following problem illustrates an unattainable solution:

\begin{displaymath}\nonumber{}\begin{array}{lclcl}\nonumber{}\mbox{minimize} & x^{2}y &  & \\\nonumber{}\mbox{subject to} & xy & \leq{} & 1,\\\nonumber{} & x,y>0. &  &\end{array}\end{displaymath}

Clearly, the optimal objective value is 0 but because of the constraint the x,y>0 constraint this value can never be attained: To see why this is a problem, remember that MOSEK substitutes [[MathCmd 832]] and [[MathCmd 833]] and solves the problem as

\begin{displaymath}\nonumber{}\begin{array}{lclcl}\nonumber{}\mbox{minimize} & e^{{2t_{x}}}e^{{t_{y}}} &  & \\\nonumber{}\mbox{subject to} & e^{{t_{x}}}e^{{t_{y}}} & \leq{} & 1,\\\nonumber{} & t_{x},t_{y}\in{}\mathbb{R}. &  &\end{array}\end{displaymath}

The optimal solution implies that [[MathCmd 835]] or [[MathCmd 836]], and thus it is unattainable.

Now, the issue should be clear: If a variable x appears only with nonnegative exponents, then fixing x=0 will minimize all terms in which it appears — but such a solution cannot be attained.

13.2.4.2. Infinite solution

A similar problem will occur if a finite optimal objective value requires a variable to be infinite. This can be illustrated by the following example:

\begin{displaymath}\nonumber{}\begin{array}{lclcl}\nonumber{}\mbox{minimize} & x^{{-2}} &  & \\\nonumber{}\mbox{subject to} & x^{{-1}} & \leq{} & 1,\\\nonumber{} & x>0, &  &\end{array}\end{displaymath}

which is a valid geometric programming problem. In this case the optimal objective is 0, but this requires x=, which is unattainable.

Again, this specific case will appear if a variable x appears only with negative exponents in the problem, implying that each term in which it appears can be minimized for [[MathCmd 838]].

13.2.5. An example

Consider the example

\begin{displaymath}\nonumber{}\begin{array}{lclcl}\nonumber{}\mbox{minimize} & x^{{-1}}y &  & \\\nonumber{}\mbox{subject to} & x^{2}y^{{-\frac{1}{2}}}+3y^{{\frac{1}{2}}}z^{{-1}} & \leq{} & 1,\\\nonumber{} & xy^{{-1}} & = & z^{2},\\\nonumber{} & -x & \leq{} & -\frac{1}{10},\\\nonumber{} & x & \leq{} & 3,\\\nonumber{} & x,y,z>0, &  &\end{array}\end{displaymath}

which is not a geometric optimization problem. However, using the obvious transformations we obtain the problem

\begin{math}\nonumber{}\begin{array}{lclcl}\nonumber{}\mbox{minimize} & x^{{-1}}y &  & \\\nonumber{}\mbox{subject to} & x^{2}y^{{-\frac{1}{2}}}+3y^{{\frac{1}{2}}}z^{{-1}} & \leq{} & 1,\\\nonumber{} & xy^{{-1}}z^{{-2}} & \leq{} & 1,\\\nonumber{} & x^{{-1}}yz^{2} & \leq{} & 1,\\\nonumber{} & \frac{1}{10}x^{{-1}} & \leq{} & 1,\\\nonumber{} & \frac{1}{3}x & \leq{} & 1,\\\nonumber{} & x,y,z>0, &  &\end{array}\end{math} (13.2.5)

which is a geometric optimization problem.

13.2.6. Solving from the command line tool

MOSEK provides the command line tool mskexpopt to solve a problem on the form (13.2.4). As demonstrated previously an optimal solution to this problem can be transformed into an optimal solution to the geometric optimization problem (13.2.1) by using the transform:

\begin{displaymath}\nonumber{}t_{j}=e^{{x_{j}}}.\end{displaymath}

A more detailed description of mskexpopt and the definition of the input format used is found in Section 6.2. The source code is also included in the MOSEK distribution.

13.2.6.1. An example

The problem (13.2.5) can be written in the mskexpopt format as follows:

5   * numcon 
3   * numvar 
7   * numter 
* Coefficients of terms
1
1
3
1
1
0.1
0.333333
* Constraints each term belong to
0
1
1
2
3
4
5
* Section defining a_kj. 
* Format: term var coef
0 0 -1
0 1  1
1 0  2
1 1  -0.5
2 1  0.5
2 2  -1
3 0 1
3 1 -1
3 2 -2
4 0 -1
4 1 1
4 2 2
5 0 -1
6 0 1 
 

The command line:

mskexpopt go1.eo

solves the problem and writes the solution file:

PROBLEM STATUS      : PRIMAL_AND_DUAL_FEASIBLE
SOLUTION STATUS     : OPTIMAL
OBJECTIVE           : 1.001904e-03

PRIMAL VARIABLES
INDEX   ACTIVITY
1       -2.302585e+00
2       -9.208438e+00
3       3.452927e+00

DUAL VARIABLES
INDEX   ACTIVITY
1       1.000000e+00
2       2.003813e+00
3       1.906415e-03
4       5.272269e+00
5       5.273223e+00
6       3.006672e+00
7       8.758884e-12

The primal solution can be transformed into a solution to the geometric optimization problem as follows

\begin{math}\nonumber{}t_{0}=e^{{-2.302585e+00}}=0.1\end{math} (13.2.6)
\begin{math}\nonumber{}t_{1}=e^{{-9.208438e+00}}=1.0019^{{-4}}\end{math} (13.2.7)
\begin{math}\nonumber{}t_{1}=e^{{3.452927e+00}}=31.5927.\end{math} (13.2.8)

13.2.7. Further information

More information about geometric optimization problems is located in [18, 24, 14].

Wed Feb 29 16:08:51 2012