7.3 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

(7.9)\[\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

(7.10)\[\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.

7.3.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\):

(7.11)\[\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 (7.11) becomes

(7.12)\[\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 (7.10). It can be written as:

(7.13)\[\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. The explicit representation of conic constraints in this case is:

\[\begin{split}\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 \EXP\times\EXP.\end{split}\]

We can now use this representation to assemble all constraints in the model.

Listing 7.8 Source code solving problem (7.12). Click here to download.

    Awall  = 200
    Afloor = 50
    alpha  = 2
    beta   = 10
    gamma  = 2
    delta  = 10
 
    % A model with variables [x,y,z,u1,u2]
    model = mosekmodel(...
        name = "gp1",...
        numvar = 5);
    model.objective("maximize", [ 1 1 1 0 0]);

    % u1 + u2 = 1
    model.appendcons(F = [ 0 0 0 1 1 ], domain = mosekdomain("equal",     rhs=1.0));
    
    % y + z <= log(Afloor)
    model.appendcons(F = [ 0 1 1 0 0 ], domain = mosekdomain("less than", rhs=log(Afloor)));
    
    % Two-sided bounds on x-y and z-y
    model.appendcons(F = [ 1 -1 0 0 0 ; ...
                           0 -1 1 0 0 ], domain = mosekdomain("greater than",dim=2, rhs=[log(alpha) log(gamma)]'));
    model.appendcons(F = [ 1 -1 0 0 0 ; ...
                           0 -1 1 0 0 ], domain = mosekdomain("less than",   dim=2, rhs=[log(beta)  log(delta)]'));

    % Conic constraints
    model.appendcons(F = sparse([1 3 3],[4 1 2],[1.0 1.0 1.0]), g = [0 1 log(alpha/Awall)]', domain = mosekdomain("exp"));
    model.appendcons(F = sparse([1 3 3],[5 1 3],[1.0 1.0 1.0]), g = [0 1 log(alpha/Awall)]', domain = mosekdomain("exp"));
    
    model.solve();
    if model.hassolution("interior")
        [xx,prosta,solsta] = model.getsolution("interior","x");
        x = xx(1:3);
        disp(exp(x));
    end