6.6 Geometric Programming

Geometric programs (GP) are a particular class of optimization problems which can be expressed in special polynomial form as positive sums of generalized monomials. More precisely, a geometric problem in canonical form is

(6.17)$$\begin{array}{r}\begin{array}{lll}\text{minimize}& f_0(x)& \\ \text{subject to}& f_i(x)\le 1,& i=1,\dots ,m,\\ & x_j>0,& j=1,\dots ,n,\end{array}\end{array}$$

where each $f_0,\dots ,f_m$ is a posynomial, that is a function of the form

$$f(x)=\sum _k{c}_k{x}_1^{{\alpha }_{k1}}{x}_2^{{\alpha }_{k2}}\cdots x_n^{{\alpha }_{kn}}$$

with arbitrary real ${\alpha }_{ki}$ and $c_k>0$. The standard way to formulate GPs in convex form is to introduce a variable substitution

$$x_i=\exp (y_i).$$

Under this substitution all constraints in a GP can be reduced to the form

(6.18)$$\log (\sum _k\exp (a_k^{T}y+b_k))\le 0$$

involving a log-sum-exp bound. Moreover, constraints involving only a single monomial in $x$ can be even more simply written as a linear inequality:

$$a_k^{T}y+b_k\le 0$$

We refer to the MOSEK Modeling Cookbook and to [BKVH07] for more details on this reformulation. A geometric problem formulated in convex form can be entered into MOSEK with the help of exponential cones.

6.6.1 Example GP1

The following problem comes from [BKVH07]. Consider maximizing the volume of a $h\times w\times d$ box subject to upper bounds on the area of the floor and of the walls and bounds on the ratios $h/w$ and $d/w$:

(6.19)$$\begin{array}{r}\begin{array}{rrl}\text{maximize}& hwd& \\ \text{subject to}& 2(hw+hd)& \le A_{\mathrm{wall}},\\ & wd& \le A_{\mathrm{floor}},\\ & \alpha & \le h/w\le \beta ,\\ & \gamma & \le d/w\le \delta .\end{array}\end{array}$$

The decision variables in the problem are $h,w,d$. We make a substitution

$$h=\exp (x),w=\exp (y),d=\exp (z)$$

after which (6.19) becomes

(6.20)$$\begin{array}{r}\begin{array}{rl}\text{maximize}& x+y+z\\ \text{subject to}& \log (\exp (x+y+\log (2/A_{\mathrm{wall}}))+\exp (x+z+\log (2/A_{\mathrm{wall}})))\le 0,\\ & y+z\le \log (A_{\mathrm{floor}}),\\ & \log (\alpha )\le x-y\le \log (\beta ),\\ & \log (\gamma )\le z-y\le \log (\delta ).\end{array}\end{array}$$

Next, we demonstrate how to implement a log-sum-exp constraint (6.18). It can be written as:

(6.21)$$\begin{array}{r}\begin{array}{l}{u}_k\ge \exp (a_k^{T}y+b_k),(\mathrm{equiv}.(u_k,1,a_k^{T}y+b_k)\in K_{\exp }),\\ \sum _k{u}_k=1.\end{array}\end{array}$$

This presentation requires one extra variable $u_k$ for each monomial appearing in the original posynomial constraint. In this case the affine conic constraints (ACC, see Sec. 6.2 (From Linear to Conic Optimization)) take the form:

$$\begin{array}{r}\left[\begin{array}{ccccc}0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0\\ 1& 1& 0& 0& 0\\ 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0\\ 1& 0& 1& 0& 0\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\\ u_1\\ u_2\end{array}\right]+\left[\begin{array}{c}0\\ 1\\ \log (2/A_{\mathrm{wall}})\\ 0\\ 1\\ \log (2/A_{\mathrm{wall}})\end{array}\right]\in K_{\exp }\times K_{\exp }.\end{array}$$

As a matter of demonstration we will also add the constraint

$$u_1+u_2-1=0$$

as an affine conic constraint. It means that to define the all the ACCs we need to produce the following affine expressions (AFE) and store them:

$$u_1,u_2,x+y+\log (2/A_{\mathrm{wall}}),x+z+\log (2/A_{\mathrm{wall}}),1.0,u_1+u_2-1.0.$$

We implement it by adding all the affine expressions (AFE) and then picking the ones required for each ACC:

Listing 6.7 Implementation of log-sum-exp as in (6.21). Click here to download.

    MSKint64t acc1_afeidx[] = {0, 4, 2};
    MSKint64t acc2_afeidx[] = {1, 4, 3};
    MSKint64t acc3_afeidx[] = {5};

    // 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
    if (r == MSK_RES_OK) 
      r = MSK_appendvars(task, 2);
    if (r == MSK_RES_OK)
      r = MSK_putvarboundsliceconst(task, numvar, numvar+2, MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);

    if (r == MSK_RES_OK) 
      r = MSK_putafefentrylist(task, f_nnz, afeidx, varidx, f_val);
    if (r == MSK_RES_OK) 
      r = MSK_putafegslice(task, 0, numafe, g);
    
    /* Append the primal exponential cone domain */
    if (r == MSK_RES_OK)
      r = MSK_appendprimalexpconedomain(task, &expdomidx);
    
    /* (u1, 1, x+y+log(2/Awall)) \in EXP */
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, expdomidx, 3, acc1_afeidx, NULL);

    /* (u2, 1, x+z+log(2/Awall)) \in EXP */
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, expdomidx, 3, acc2_afeidx, NULL);

    /* The constraint u1+u2-1 \in \ZERO is added also as an ACC */
    if (r == MSK_RES_OK)
      r = MSK_appendrzerodomain(task, 1, &rzerodomidx);
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, rzerodomidx, 1, acc3_afeidx, NULL);

We can now use this function to assemble all constraints in the model. The linear part of the problem is entered as in Sec. 6.1 (Linear Optimization).

Listing 6.8 Source code solving problem (6.20). Click here to download.

int max_volume_box(double Aw, double Af, 
                   double alpha, double beta, double gamma, double delta,
                   double hwd[])
{
  // Basic dimensions of our problem
  const int numvar    = 3;  // Variables in original problem
  const int numcon    = 3;  // Linear constraints in original problem

  // Linear part of the problem involving x, y, z
  const double       cval[]  = {1, 1, 1};
  const int          asubi[] = {0, 0, 1, 1, 2, 2};
  const int          asubj[] = {1, 2, 0, 1, 2, 1};
  const int          alen    = 6;
  const double       aval[]  = {1.0, 1.0, 1.0, -1.0, 1.0, -1.0};
  const MSKboundkeye bkc[]   = {MSK_BK_UP, MSK_BK_RA, MSK_BK_RA};
  const double       blc[]   = {-MSK_INFINITY, log(alpha), log(gamma)};
  const double       buc[]   = {log(Af), log(beta), log(delta)};

  // Affine conic constraint data of the problem
  MSKint64t       expdomidx, rzerodomidx;
  const MSKint64t numafe = 6, f_nnz = 8;
  const MSKint64t afeidx[] = {0, 1, 2, 2, 3, 3, 5, 5};
  const MSKint32t varidx[] = {3, 4, 0, 1, 0, 2, 3, 4};
  const double     f_val[] = {1, 1, 1, 1, 1, 1, 1, 1};
  const double         g[] = {0, 0, log(2/Aw), log(2/Aw), 1, -1};

  MSKtask_t          task = NULL;
  MSKrescodee        r = MSK_RES_OK, trmcode;
  MSKsolstae         solsta;
  MSKint32t          i;
  double             *xyz = (double*) calloc(numvar, sizeof(double));

  if (r == MSK_RES_OK)
    r = MSK_maketask(NULL, 0, 0, &task);

  if (r == MSK_RES_OK)
    r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

  if (r == MSK_RES_OK)
    r = MSK_appendvars(task, numvar);

  if (r == MSK_RES_OK)
    r = MSK_appendcons(task, numcon);

  if (r == MSK_RES_OK)
    r = MSK_appendafes(task, numafe);

  // Objective is the sum of three first variables
  if (r == MSK_RES_OK)
    r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);
  if (r == MSK_RES_OK)
    r = MSK_putcslice(task, 0, numvar, cval);
  if (r == MSK_RES_OK)
    r = MSK_putvarboundsliceconst(task, 0, numvar, MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);

  // Add the three linear constraints
  if (r == MSK_RES_OK)
    r = MSK_putaijlist(task, alen, asubi, asubj, aval);
  if (r == MSK_RES_OK)  
    r = MSK_putconboundslice(task, 0, numvar, bkc, blc, buc);

  if (r == MSK_RES_OK)
  {
    MSKint64t acc1_afeidx[] = {0, 4, 2};
    MSKint64t acc2_afeidx[] = {1, 4, 3};
    MSKint64t acc3_afeidx[] = {5};

    // 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
    if (r == MSK_RES_OK) 
      r = MSK_appendvars(task, 2);
    if (r == MSK_RES_OK)
      r = MSK_putvarboundsliceconst(task, numvar, numvar+2, MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);

    if (r == MSK_RES_OK) 
      r = MSK_putafefentrylist(task, f_nnz, afeidx, varidx, f_val);
    if (r == MSK_RES_OK) 
      r = MSK_putafegslice(task, 0, numafe, g);
    
    /* Append the primal exponential cone domain */
    if (r == MSK_RES_OK)
      r = MSK_appendprimalexpconedomain(task, &expdomidx);
    
    /* (u1, 1, x+y+log(2/Awall)) \in EXP */
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, expdomidx, 3, acc1_afeidx, NULL);

    /* (u2, 1, x+z+log(2/Awall)) \in EXP */
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, expdomidx, 3, acc2_afeidx, NULL);

    /* The constraint u1+u2-1 \in \ZERO is added also as an ACC */
    if (r == MSK_RES_OK)
      r = MSK_appendrzerodomain(task, 1, &rzerodomidx);
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, rzerodomidx, 1, acc3_afeidx, NULL);
  }

  // Solve and map to original h, w, d
  if (r == MSK_RES_OK)
    r = MSK_optimizetrm(task, &trmcode);

  if (r == MSK_RES_OK)
    MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

  if (solsta == MSK_SOL_STA_OPTIMAL)
  {
    if (r == MSK_RES_OK)
      r = MSK_getxxslice(task, MSK_SOL_ITR, 0, numvar, xyz);
    for(i = 0; i < numvar; i++) hwd[i] = exp(xyz[i]);
  }
  else
  {
    printf("Solution not optimal, termination code %d.\n", trmcode);
    r = trmcode;
  }

  free(xyz);
  return r;
}

Given sample data we obtain the solution $h,w,d$ as follows:

Listing 6.9 Sample data for problem (6.19). Click here to download.

int main()
{
  const double Aw    = 200.0;
  const double Af    = 50.0;
  const double alpha = 2.0;
  const double beta  = 10.0;
  const double gamma = 2.0;
  const double delta = 10.0;
  MSKrescodee  r;
  double       hwd[3];

  r = max_volume_box(Aw, Af, alpha, beta, gamma, delta, hwd);

  printf("Response code: %d\n", r);
  if (r == MSK_RES_OK)
    printf("Solution h=%.4f w=%.4f d=%.4f\n", hwd[0], hwd[1], hwd[2]);

  return r;
}