6.8 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{split}\begin{array}{lll} \minimize & f_0(x) & \\ \st & f_i(x) \leq 1, &i=1,\ldots,m, \\ & x_j>0, &j=1,\ldots,n, \end{array}\end{split}\]

where each \(f_0,\ldots,f_m\) is a posynomial, that is a function of the form

\[f(x) = \sum_k c_kx_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^Ty+b_k)) \leq 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^Ty+b_k\leq 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.8.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{split}\begin{array}{rrl} \maximize & hwd & \\ \st & 2(hw + hd) & \leq A_{\mathrm{wall}}, \\ & wd & \leq A_{\mathrm{floor}}, \\ & \alpha & \leq h/w \leq \beta, \\ & \gamma & \leq d/w \leq \delta. \end{array}\end{split}\]

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{split}\begin{array}{rll} \maximize & x+y+z \\ \st & \log(\exp(x+y+\log(2/A_{\mathrm{wall}}))+\exp(x+z+\log(2/A_{\mathrm{wall}}))) \leq 0, \\ & y+z \leq \log(A_{\mathrm{floor}}), \\ & \log(\alpha) \leq x-y \leq \log(\beta), \\ & \log(\gamma) \leq z-y \leq \log(\delta). \end{array}\end{split}\]

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

(6.21)\[\begin{split}\begin{array}{l} u_k\geq \exp(a_k^Ty+b_k),\quad (\mathrm{equiv.}\ (u_k,1,a_k^Ty+b_k)\in\EXP),\\ \sum_k u_k = 1. \end{array}\end{split}\]

This presentation requires one extra variable \(u_k\) for each monomial appearing in the original posynomial constraint. Another fixed variable \(t_k=1\) stands for the second entry in the exponential cone.

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

          // Add a single log-sum-exp constraint sum(log(exp(z_i))) <= 0
          // Assume numExp variable triples are ordered as (u0,t0,z0,u1,t1,z1...)
          // starting from variable with index expStart
          double[]  val = new double[numExp];
          int[]     sub = new int[numExp];

          // sum(u_i) = 1 as constraint number c, u_i unbounded
          for(int i = 0; i < numExp; i++)
            { sub[i] = expStart + 3*i; val[i] = 1.0; }
          task.putarow(c, sub, val);
          task.putconbound(c, boundkey.fx, 1.0, 1.0);
          task.putvarboundlistconst(sub, boundkey.fr, -inf, inf);

          // z_i unbounded
          for(int i = 0; i < numExp; i++) sub[i] = expStart + 3*i + 2;
          task.putvarboundlistconst(sub, boundkey.fr, -inf, inf);

          // t_i = 1
          for(int i = 0; i < numExp; i++) sub[i] = expStart + 3*i + 1; 
          task.putvarboundlistconst(sub, boundkey.fx, 1, 1);

          // Every triple is in an exponential cone
          conetype[] ct = new conetype[numExp];
          int[]     len = new int[numExp];
          for(int i = 0; i < numExp; i++) 
            { ct[i] = conetype.pexp; val[i] = 0.0; len[i] = 3; }
          task.appendconesseq(ct, val, len, expStart);

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.14 Source code solving problem (6.20). Click here to download.

    public static double[] max_volume_box(double Aw, double Af, 
                                          double alpha, double beta, double gamma, double delta)
    {
      // Basic dimensions of our problem
      int numvar    = 3;  // Variables in original problem
      int numLinCon = 3;  // Linear constraints in original problem
      int numExp    = 2;  // Number of exp-terms in the log-sum-exp constraint

      // Linear part of the problem
      double[] cval  = {1, 1, 1};
      int[]    asubi = {0, 0, 1, 1, 2, 2};
      int[]    asubj = {1, 2, 0, 1, 2, 1};
      double[] aval  = {1.0, 1.0, 1.0, -1.0, 1.0, -1.0};
      boundkey[] bkc = {boundkey.up, boundkey.ra, boundkey.ra};
      double[] blc   = {-inf, Math.Log(alpha), Math.Log(gamma)};
      double[] buc   = {Math.Log(Af), Math.Log(beta), Math.Log(delta)};

      // Linear part setting up slack variables
      // for the linear expressions appearing inside exps
      // x_5 - x - y = log(2/Awall)
      // x_8 - x - z = log(2/Awall)
      // The slack indexes are convenient for defining exponential cones, see later
      int[]    a2subi = {3, 3, 3, 4, 4, 4};
      int[]    a2subj = {5, 0, 1, 8, 0, 2};
      double[] a2val  = {1.0, -1.0, -1.0, 1.0, -1.0, -1.0};
      boundkey[] b2kc = {boundkey.fx, boundkey.fx};
      double[] b2luc  = {Math.Log(2/Aw), Math.Log(2/Aw)};

      using (Env env = new Env())
      {
        using (Task task = new Task(env, 0, 0))
        {
          // Directs the log task stream to the user specified
          // method task_msg_obj.stream
          task.set_Stream (mosek.streamtype.log, new msgclass (""));

          // Add variables and constraints
          task.appendvars(numvar + 3*numExp);
          task.appendcons(numLinCon + numExp + 1);

          // 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);

          // Add linear constraints for the expressions appearing in exp(...)
          task.putaijlist(a2subi, a2subj, a2val);
          task.putconboundslice(numLinCon, numLinCon+numExp, b2kc, b2luc, b2luc);

          int c = numLinCon + numExp;
          int expStart = numvar;
          // Add a single log-sum-exp constraint sum(log(exp(z_i))) <= 0
          // Assume numExp variable triples are ordered as (u0,t0,z0,u1,t1,z1...)
          // starting from variable with index expStart
          double[]  val = new double[numExp];
          int[]     sub = new int[numExp];

          // sum(u_i) = 1 as constraint number c, u_i unbounded
          for(int i = 0; i < numExp; i++)
            { sub[i] = expStart + 3*i; val[i] = 1.0; }
          task.putarow(c, sub, val);
          task.putconbound(c, boundkey.fx, 1.0, 1.0);
          task.putvarboundlistconst(sub, boundkey.fr, -inf, inf);

          // z_i unbounded
          for(int i = 0; i < numExp; i++) sub[i] = expStart + 3*i + 2;
          task.putvarboundlistconst(sub, boundkey.fr, -inf, inf);

          // t_i = 1
          for(int i = 0; i < numExp; i++) sub[i] = expStart + 3*i + 1; 
          task.putvarboundlistconst(sub, boundkey.fx, 1, 1);

          // Every triple is in an exponential cone
          conetype[] ct = new conetype[numExp];
          int[]     len = new int[numExp];
          for(int i = 0; i < numExp; i++) 
            { ct[i] = conetype.pexp; val[i] = 0.0; len[i] = 3; }
          task.appendconesseq(ct, val, len, expStart);

          // Solve and map to original h, w, d
          task.optimize();
          double[] xyz = new double[numvar];
          double[] hwd = new double[numvar];
          task.getxxslice(soltype.itr, 0, numvar, xyz);
          for(int i = 0; i < numvar; i++) hwd[i] = Math.Exp(xyz[i]);
          return hwd;
        }
      }
    }

Given sample data we obtain the solution \(h,w,d\) as follows:

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

    public static void Main(String[] args)
    {
      double Aw    = 200.0;
      double Af    = 50.0;
      double alpha = 2.0;
      double beta  = 10.0;
      double gamma = 2.0;
      double delta = 10.0;
      
      double[] hwd = max_volume_box(Aw, Af, alpha, beta, gamma, delta);

      Console.WriteLine("h={0:f4} w={1:f4} d={2:f4}", hwd[0], hwd[1], hwd[2]);
    }