Geometric programming

A geometric program (GP) is a constrained optimization problem built from monomials and posynomials. A monomial is a function

\[ g(x) = c\,x_1^{a_1} x_2^{a_2} \cdots x_n^{a_n}, \qquad c > 0, \;\; a_i \in \mathbb{R}, \]

with a strictly positive coefficient and real exponents (so \(x/y\) and \(\sqrt{x}\) are monomials, but \(x - y\) is not). A posynomial is a sum of monomials. A GP in standard form minimizes a posynomial subject to posynomial upper bounds and monomial equalities [Boyd et al., 2007]:

\[\begin{split} \begin{aligned} \text{minimize} \quad & f_0(x) \\ \text{subject to} \quad & f_i(x) \le 1, && i = 1,\dots,m \quad (\text{posynomial}) \\ & g_j(x) = 1, && j = 1,\dots,p \quad (\text{monomial}) \\ & x > 0. \end{aligned} \end{split}\]

A posynomial is not convex in \(x\) — its Hessian is indefinite on the positive orthant. The crucial fact, going back to Duffin et al. [1967], is that the logarithmic change of variables \(u_i = \log x_i\) makes the problem convex. Under \(x_i = e^{u_i}\) a monomial \(c\prod_i x_i^{a_i}\) becomes \(\exp(a^\top u + \log c)\), a posynomial becomes a sum of exponentials of affine forms (a convex function), the constraint \(f_i(x)\le 1\) becomes a convex sublevel set, and a monomial equality becomes a plain linear equality. So a GP is a convex program in disguise: it has a unique global optimum reachable by a single convex NLP solve, with no branch-and-bound and no spurious local minima [Boyd and Vandenberghe, 2004].

This log-convexity is why GPs are a workhorse of engineering design — analog circuit sizing, structural and mechanical design, and aspect-ratio / packing problems all fit the monomial/posynomial template [Boyd et al., 2007]. discopt’s discopt.gp module recognizes GP structure in an ordinary Model, builds the equivalent convex log-space formulation, and solves it. The same machinery that powers disciplined geometric programming in modeling frameworks like CVXPY [Agrawal et al., 2019] underlies this recognition.

The three entry points are:

• classify_gp(model) — returns a GPStructure if the model is a GP, else None.

• as_geometric_program(model) — builds the convex log-space NLP plus the \(x = \exp(u)\) recovery map (or None).

• solve_gp(model, **kwargs) — solves via the log-space transform and maps the solution back to \(x\)-space (or None for a non-GP).

import discopt.modeling as dm
from discopt.gp import as_geometric_program, classify_gp, solve_gp

A first geometric program

We build a tiny GP in two variables. Note the two GP-specific requirements: every variable must be continuous and strictly positive (we give a small positive lower bound), and every term must be a monomial with a positive coefficient. The objective \(x/y + 2xy\) is a posynomial, and \(xy + x/y \le 10\) is a posynomial bounded above by a constant — both legal GP ingredients.

We also add \(1/x \le 1\) (a monomial \(\le 1\), i.e. \(x \ge 1\)), and it is worth saying why. Write the objective as \(x\,(1/y + 2y)\): for any fixed \(y\) it is strictly increasing in \(x\), so without a real lower bound on \(x\) the infimum is \(0\) and the solver would simply drive \(x\) onto whatever lb we typed. The lb=1e-4 on the variable is a numerical device to keep the log-transform well defined, not a modeling statement — an optimum sitting on it is an artifact of the device, not an answer. \(1/x \le 1\) is the genuine bound.

With it the problem has a closed-form solution we can check: minimizing \(1/y + 2y\) gives \(y^\star = 1/\sqrt{2}\), and \(x^\star = 1\) at the bound, so \(f^\star = 2\sqrt{2} \approx 2.8284\).

m = dm.Model("gp1")
x = m.continuous("x", lb=1e-4, ub=1e4)
y = m.continuous("y", lb=1e-4, ub=1e4)

m.minimize(x / y + 2.0 * x * y)
m.subject_to(x * y + x / y <= 10.0)  # posynomial <= constant : a valid GP constraint
m.subject_to(1.0 / x <= 1.0)  # monomial <= 1, i.e. x >= 1 : the real lower bound

structure = classify_gp(m)
print("recognized as a GP:", structure is not None)
print("minimize         :", structure.minimize)
print("objective terms  :", len(structure.objective.monomials), "monomials")
print("GP constraints   :", len(structure.constraints))

recognized as a GP: True
minimize         : True
objective terms  : 2 monomials
GP constraints   : 2

classify_gp returned a GPStructure, so discopt has a proof that this model is log-convex. We can inspect the convex reformulation it builds: a new model in the log-variables \(u = \log x\), where each posynomial becomes a sum of exponentials of affine forms.

gp = as_geometric_program(m)
print("log-space scalar variables:", gp.n_scalars)
print("log-space model name      :", gp.log_model.name)
print("log-space constraints     :", len(gp.log_model._constraints))

log-space scalar variables: 2
log-space model name      : gp1_log
log-space constraints     : 2

Solving it

solve_gp solves the convex log-space NLP and maps the optimum back to the original \(x\)-space. Because the transformed program is convex and equivalent to the GP, an optimal status certifies the global optimum — the reported objective is simultaneously the best incumbent and a valid lower bound, so the optimality gap is exactly zero with no branch-and-bound.

import math

result = solve_gp(m)
print(f"status   : {result.status}")
print(f"objective: {float(result.objective):.6f}")
print(f"gap      : {result.gap}")
for name, value in result.x.items():
    print(f"  {name} = {float(value):.6f}")

# Check against the closed form derived above, rather than eyeballing it.
print(f"\nanalytic : x*=1, y*=1/sqrt(2)={1 / math.sqrt(2):.6f}, "
      f"f*=2*sqrt(2)={2 * math.sqrt(2):.6f}")
assert abs(float(result.objective) - 2 * math.sqrt(2)) < 1e-5, result.objective
assert abs(float(result.x["y"]) - 1 / math.sqrt(2)) < 1e-5, result.x["y"]

# And confirm the answer is NOT an artifact of the variable's numerical lb:
# re-solving with a far smaller lb must give the same optimum.
m2 = dm.Model("gp1_tighter_lb")
x2 = m2.continuous("x", lb=1e-8, ub=1e4)
y2 = m2.continuous("y", lb=1e-8, ub=1e4)
m2.minimize(x2 / y2 + 2.0 * x2 * y2)
m2.subject_to(x2 * y2 + x2 / y2 <= 10.0)
m2.subject_to(1.0 / x2 <= 1.0)
print(f"same model with lb=1e-8: objective = {float(solve_gp(m2).objective):.6f}")

status   : optimal
objective: 2.828427
gap      : 0.0
  x = 1.000000
  y = 0.707107

analytic : x*=1, y*=1/sqrt(2)=0.707107, f*=2*sqrt(2)=2.828427
same model with lb=1e-8: objective = 2.828428

A worked design example: maximum-volume box

A classic GP from the engineering-design literature is sizing a box to maximize its volume subject to limits on the total wall and floor area and on its aspect ratios [Boyd et al., 2007]. With height \(h\), width \(w\), and depth \(d\),

\[\begin{split} \begin{aligned} \text{maximize} \quad & h\,w\,d \\ \text{s.t.} \quad & 2(h w + h d) \le A_\text{wall}, \\ & w d \le A_\text{floor}, \\ & \alpha \le h/w \le \beta, \qquad \gamma \le d/w \le \delta. \end{aligned} \end{split}\]

Maximizing \(hwd\) is the same as minimizing the monomial \(1/(hwd)\), which keeps the objective a monomial. Each area limit divides through by its constant to give a posynomial \(\le 1\), and each two-sided aspect-ratio bound splits into two monomial \(\le 1\) constraints. The result is a genuine GP. (discopt parses each summand of a posynomial individually, so we write the wall-area constraint as an explicit sum of monomials rather than a scalar times a parenthesized sum.)

A_wall, A_floor = 100.0, 10.0
alpha, beta = 0.5, 2.0  # 0.5 <= h/w <= 2.0
gamma, delta = 0.5, 2.0  # 0.5 <= d/w <= 2.0

box = dm.Model("box")
h = box.continuous("h", lb=1e-3, ub=1e3)
w = box.continuous("w", lb=1e-3, ub=1e3)
d = box.continuous("d", lb=1e-3, ub=1e3)

box.minimize(1.0 / (h * w * d))  # maximize volume  <=>  minimize 1/volume

box.subject_to((2.0 / A_wall) * h * w + (2.0 / A_wall) * h * d <= 1.0)  # wall area
box.subject_to((1.0 / A_floor) * w * d <= 1.0)  # floor area
box.subject_to(alpha * w / h <= 1.0)  # alpha <= h/w
box.subject_to((1.0 / beta) * h / w <= 1.0)  # h/w <= beta
box.subject_to(gamma * w / d <= 1.0)  # gamma <= d/w
box.subject_to((1.0 / delta) * d / w <= 1.0)  # d/w <= delta

print("recognized as a GP:", classify_gp(box) is not None)

recognized as a GP: True

res = solve_gp(box)
hv, wv, dv = (float(res.x[name]) for name in ("h", "w", "d"))

print(f"status      : {res.status}")
print(f"dimensions  : h={hv:.4f}  w={wv:.4f}  d={dv:.4f}")
print(f"volume      : {hv * wv * dv:.4f}")
print(f"wall area   : {2 * (hv * wv + hv * dv):.4f}  (limit {A_wall})")
print(f"floor area  : {wv * dv:.4f}  (limit {A_floor})")
print(f"aspect h/w  : {hv / wv:.4f}   d/w : {dv / wv:.4f}")

status      : optimal
dimensions  : h=7.7460  w=3.8730  d=2.5820
volume      : 77.4597
wall area   : 100.0000  (limit 100.0)
floor area  : 10.0000  (limit 10.0)
aspect h/w  : 2.0000   d/w : 0.6667

The optimum drives the wall- and floor-area constraints to their limits and lands on the \(h/w\) aspect bound — exactly the behavior we expect from a well-posed design GP. Because the log-space problem is convex, this solution is the certified global maximum-volume box.

Monomial equalities are allowed

GP standard form permits monomial equality constraints (a monomial equal to a constant, or to another monomial). The constraint \(xy = 4\) couples the two variables along a hyperbola; under \(u = \log x\) it becomes the linear equality \(u_x + u_y = \log 4\), which keeps the log-space problem convex.

eq = dm.Model("gp_eq")
a = eq.continuous("a", lb=1e-4, ub=1e4)
b = eq.continuous("b", lb=1e-4, ub=1e4)

eq.minimize(a + b)
eq.subject_to(a * b == 4.0)  # monomial == constant : a valid GP equality

print("recognized as a GP:", classify_gp(eq) is not None)
eq_res = solve_gp(eq)
print(f"status   : {eq_res.status}")
print(f"objective: {float(eq_res.objective):.6f}   (a+b minimized at a=b=2)")
print(f"a={float(eq_res.x['a']):.4f}  b={float(eq_res.x['b']):.4f}")

recognized as a GP: True
status   : optimal
objective: 4.000000   (a+b minimized at a=b=2)
a=2.0000  b=2.0000

A non-example: signomials

GP recognition is conservative: discopt only declares a model a GP when it can prove every term is a monomial with a positive coefficient. A term with a negative coefficient breaks the posynomial structure — the result is a signomial, for which the log-transform is no longer convex [Boyd et al., 2007].

The constraint \(x y - 3x \le 1\) contains the negative-coefficient term \(-3x\). There is no way to write it as posynomial \(\le\) monomial, so classify_gp returns None and solve_gp declines the model rather than silently solving a problem it cannot certify.

sig = dm.Model("signomial")
p = sig.continuous("p", lb=1e-4, ub=1e4)
q = sig.continuous("q", lb=1e-4, ub=1e4)

sig.minimize(p + q)
sig.subject_to(p * q - 3.0 * p <= 1.0)  # negative coefficient -> signomial, NOT a GP

print("classify_gp returns None:", classify_gp(sig) is None)
print("solve_gp returns None   :", solve_gp(sig) is None)

classify_gp returns None: True
solve_gp returns None   : True

When to reach for GP detection

GP recognition is an additive fast path: when a model happens to be a geometric program, discopt can solve it to global optimality with a single convex NLP solve, sidestepping spatial branch-and-bound entirely. This is a large win for the design problems that GPs were built for — circuit sizing, structural and mechanical proportions, and packing/aspect-ratio problems [Boyd et al., 2007].

A few practical notes:

• Strict positivity is mandatory. Every variable needs a strictly positive lower bound, since the optimum lives in log-space. Use a small bound like lb=1e-4 rather than 0.

• Keep coefficients positive. Any term with a negative coefficient pushes the model out of GP form into the signomial class, where global optimality is no longer free. Detection is deliberately conservative so it never mis-certifies such a model.

• Detection is sound but not complete. classify_gp recognizes a model that is already written in (expanded) posynomial/monomial form; an algebraically equivalent GP written differently may not be recognized. When in doubt, call classify_gp and fall back to discopt’s general MINLP solver if it returns None.

The GP log-space verdict is kept deliberately separate from discopt’s \(x\)-space convexity detection: a genuine GP is generally not convex in \(x\), so folding log-convexity into the ordinary convex fast path would mis-gate it. The two predicates coexist, each guarding its own equivalent convex reformulation [Boyd and Vandenberghe, 2004].