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{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

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

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

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

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

$$\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}$$

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