Engineering Design by Geometric Programming

Abstract

A geometric program (GP) is a type of mathematical optimization problem characterized by objective and constraint functions, where all functions are of signomial form. The importance of GP comes from two relatively recent developments: (i) new methods can solve even large-scale GP extremely efficiently and reliably; (ii) a number of practical problems have recently been found to be equivalent to or approximated by GP. This study proposes an optimization approach for solving GP. Our approach is first to convert all signomial terms in GP into convex and concave terms. Then the concave terms are further treated with the proposed piecewise linearization method where only ⌈log₂(𝑚−1)⌉ binary variables are used. It has the following features: (i) it offers more convenient and efficient means of expressing a piecewise linear function; (ii) fewer 0-1 variables are used; (iii) the computational results show that the proposed method is much more efficient and faster than the conventional one, especially when the number of break points becomes large. In addition, the engineering design problems are illustrated to evaluate the usefulness of the proposed methods.

1. Introduction

Consider the following geometric program (GP) [1]: min𝑓₀(𝑥)s.t.𝑓_𝑘(𝑥)𝑔≤1,𝑘=1,2,…,𝑝,_𝑙(𝑥)𝑥=1,𝑙=1,2,…,𝑞,𝑥=₁,𝑥₂,…,𝑥_𝑛, (1)(2)(3)(4) where 𝑔_𝑙(𝑥), 𝑙=1,2,…,𝑞, are monomials as in (5), 𝑐_𝑙≥0, and the exponential constants 𝛼_𝑖𝑙∈ℝ, 𝑖=1,2,…,𝑛, 𝑙=1,2,…,𝑞: 𝑔_𝑙(𝑥)=𝑐_𝑙^𝑛_𝑖=1𝑥^𝛼𝑖𝑙_𝑖 (5) and 𝑓_𝑘(𝑥), 𝑘=0,1,…,𝑝, are posynomials as in (6) which are sum of one or more monomials, 𝑐_ℎ𝑘≥0, and 𝛼_𝑖ℎ𝑘∈ℝ,𝑖=1,2,…,𝑛, ℎ=1,2,…,𝐻_𝑘,𝑘=0,1,…,𝑝: 𝑓_𝑘(𝑥)=^𝐻𝑘_ℎ=1𝑐^𝑛_{ℎ𝑘𝑖=1}𝑥^{𝛼𝑖ℎ𝑘}_𝑖. (6)

GP has been applied in many fields of applications including analog/digital circuit design [2–8], chemical engineering [9–11], mechanical engineering [12–19], power control  [20], and communication network systems [21–23].

Obtaining the optimal solutions for GP is not straightforward because the signomial terms in the objective function and constraints cannot be solved directly. As a result, many approaches have been developed. Coello and Cortés [24] proposed a genetic algorithm with an artificial immune system to solve a GP in engineering optimization. Nevertheless, this method can only obtain the local optima. Horst and Tuy [25] introduced an analytical approach for solving a problem with Lipschitzian objective and constraints. The restriction of this approach is to find the global optimum only if the range of variables can be reduced by analytical techniques. Sherali and Tuncbilek [26] developed a reformulation-linearization technique (RLT) which generates polynomial implied constraints and then linearizes the resulting problem by introducing new variables. Lin and Tsai [14] introduced a generalized method to find multiple optimal solutions of signomial discrete programming problems with free variables. By means of variable substitution and convexification strategies, a signomial discrete programming problem with free variables is first converted into another convex mixed-integer nonlinear programming problem solvable to obtain an exactly global optimum. Tsai [17] proposed a novel method to solve signomial discrete programming problems. The signomial terms are first convexified following the convexification strategies. The original program is then converted into a convex integer program solvable by commercialized packages to obtain globally optimal solutions. Tsai and Lin [19] developed an efficient method to solve a posynomial geometric program with separable functions. Power transformations and exponential transformations are utilized to convexify and underestimate posynomial terms. The original program therefore can be converted into a convex mixed-integer nonlinear program solvable to obtain a global optimum.

This paper develops an optimization approach for solving GP. Our approach is listed as follows.(i)Convert all signomial terms into convex and concave ones. (ii)The concave terms are further treated with the proposed piecewise linearization method where only ⌈log₂(𝑚−1)⌉ binary variables are used.

The rest of this paper is organized as follows. Section 2 introduces the proposed methods. Section 3 provides some numerical examples to illustrate the modeling idea and the usefulness of the proposed method. Section 4 gives our conclusions.

2. Proposed Methods

In this section, the proposed methods will be presented. First, we recall what a signomial function is. A function 𝑓∶ℝ^𝑛₊₊→ℝ defined as 𝑓(𝑥)=𝑐^𝑛_𝑖=1𝑥^𝛼𝑖_𝑖, (7) where 𝑐>0 and 𝛼_𝑖∈ℝ for all 𝑖, is called a monomial function or simply a monomial [1]. Note that the exponents 𝛼_𝑖 of a monomial can be any real numbers, but the coefficient 𝑐 must be nonnegative. A sum of monomials, namely, a function of the form 𝑓(𝑥)=^𝑁_𝑘=1𝑐_𝑘^𝑛_𝑖=1𝑥^𝛼𝑖𝑘_𝑖, (8) where 𝑐_𝑘>0 and 𝛼_𝑖𝑘∈ℝ, is called a posynomial function with 𝑁 terms or simply a posynomial. A signomial is a linear combination of monomials of some positive variables 𝑥₁,…,𝑥_𝑛.

Let 𝑓 be defined on an open convex set 𝐷⊆ℝ^𝑛 and be twice differentiable; it is known that (i) 𝑓 is convex on 𝐷 if and only if the Hessian matrix ∇²𝑓(𝑥) is positive semidefinite (p.s.d. for short) at each 𝑥∈𝐷; (ii) if ∇²𝑓(𝑥) is positive definite (p.d. for short) at each 𝑥∈𝐷, then 𝑓 is strictly convex. For convenience, we denote by Δ_𝑘 the leading principal minors of 𝐴.

Lemma 1. Let 𝐴 be a nonzero 𝑛×𝑛 symmetric matrix.(a)If 𝐴 is positive semidefinite, then Δ_𝑘≥0,𝑘=1,2,…,𝑛, and not all Δ_𝑘=0.(b)𝐴 is positive definite if and only if Δ_𝑘>0, for all 𝑘=1,2,…,𝑛.

The converse of Lemma 1(a) is false. For example, let 𝐴=¹⁰⁰⁰⁰⁰₀₀₋₁; we have ⟨𝑥,𝐴𝑥⟩=𝑥²₁−𝑥²₃ which is not always nonnegative for all 𝑥∈ℝ³. But Δ₁=1≥0, Δ₂=0≥0, and Δ₃=0≥0. In fact, the converse of Lemma 1(a) is true only for 𝑛=2; see [27, page 112].

Lemma 2. Let 𝐴 be a nonzero symmetric matrix. Then, the following will hold.(a)𝐴 is 𝑝.𝑠.𝑑. if and only if all of its eigenvalues are nonnegative.(b)𝐴 is 𝑝.𝑑. if and only if all of its eigenvalues are positive.

Proposition 3. Let 𝑓₁∶ℝ^𝑛₊₊→ℝ be defined as 𝑓₁(𝑥)=𝑐₁∏^𝑛_𝑖=1𝑥^𝛼𝑖_𝑖, where 𝑐₁>0 and 𝛼_𝑖≤0 for all 𝑖=1,2,…,𝑛. Then 𝑓₁ is a convex function.

Proof. Since 𝑐₁>0, it is enough to show that 𝑓₁(𝑥)=∏^𝑛_𝑖=1𝑥^𝛼𝑖_𝑖 is convex. Let 𝑔(𝑥)=ln𝑓₁(𝑥)=∑^𝑛_𝑖=1ln𝑥^𝛼𝑖_𝑖=∑^𝑛_𝑖=1𝛼_𝑖ln𝑥_𝑖. Then, we have (𝑥)=𝛼∇𝑔¹𝑥₁𝛼₂𝑥₂⋯𝛼_𝑛𝑥_𝑛^𝑇,∇²𝑔(𝑥)=−𝛼₁𝑥²₁00⋯0−𝛼₂𝑥²₂⋯0⋮⋮⋱⋮00⋯−𝛼_𝑛𝑥²_𝑛. (9) Because 𝛼_𝑖≤0 for all 𝑖=1,2,…,𝑛, we know that all eigenvalues of ∇²𝑔(𝑥) are nonnegative, which implies (by Lemma 1(a)) that ∇²𝑔(𝑥) is positive semidefinite. Thus, 𝑔(𝑥)=ln𝑓(𝑥) is a convex function which yields that 𝑓₁(𝑥) is a convex function.

Proposition 4. Let 𝑓₂∶ℝ^𝑛₊₊→ℝ be defined as 𝑓₂(𝑥)=𝑐₂∏^𝑛_𝑖=1𝑥^𝛼𝑖_𝑖, where 𝑐₂<0 and 𝛼_𝑖>0 for all 𝑖=1,2,…,𝑛 with 1−∑^𝑛_𝑖=1𝛼_𝑖≥0. Then 𝑓₂ is a convex function.

Proof. It is not hard to compute that [∇𝑓₂(𝑥)]_𝑖=𝑐₂𝛼_𝑖𝑥^𝛼𝑖_𝑖⁻¹∏^𝑛_{𝑗=1,𝑗≠𝑖}𝑥^𝛼𝑗_𝑗. In other words, ∇𝑓₂(𝑥)=𝑐₂𝛼₁𝑥^𝛼1₁⁻¹𝑥^𝛼2₂⋯𝑥^𝛼𝑛_𝑛𝑐₂𝛼₂𝑥^𝛼1₁𝑥^𝛼2₂⁻¹⋯𝑥^𝛼𝑛_𝑛⋮𝑐₂𝛼_𝑛𝑥^𝛼1₁𝑥^𝛼2₂⋯𝑥^𝛼𝑛_𝑛⁻¹. (10) In addition, it can be verified that ∇²𝑓₂(𝑥)_𝑖𝑗=𝜕²𝑓₂(𝑥)𝜕𝑥_𝑖𝜕𝑥_𝑗=𝛼_𝑖𝛼_𝑗𝜕𝑥_𝑖𝜕𝑥_𝑗𝑓₂(𝑥)𝛼,𝑖≠𝑗,_𝑖𝛼_𝑖−1𝜕𝑥²_𝑖𝑓₂(𝑥),𝑖=𝑗, (11) Namely,∇²𝑓₂(𝑥)=𝑐₂𝛼₁𝛼₁−1𝑥₁^𝑛−2∏_𝑖=1𝑥^𝛼𝑖_𝑖𝑐₂𝛼₁𝛼₂𝑥₁⁻¹𝑥₂^𝑛−1∏_𝑖=1𝑥^𝛼𝑖_𝑖⋯𝑐₂𝛼₁𝛼_𝑛𝑥₁⁻¹𝑥_𝑛^𝑛−1∏_𝑖=1𝑥^𝛼𝑖_𝑖𝑐₂𝛼₂𝛼₁𝑥₂⁻¹𝑥₁^𝑛−1∏_𝑖=1𝑥^𝛼𝑖_𝑖𝑐₂𝛼₂𝛼₂−1𝑥₂^𝑛−2∏_𝑖=1𝑥^𝛼𝑖_𝑖⋯𝑐₂𝛼₂𝛼_𝑛𝑥₂⁻¹𝑥_𝑛^𝑛−1∏_𝑖=1𝑥^𝛼𝑖_𝑖𝑐⋮⋮⋱⋮₂𝛼_𝑛𝛼₁𝑥_𝑛⁻¹𝑥₁^𝑛−1∏_𝑖=1𝑥^𝛼𝑖_𝑖𝑐₂𝛼_𝑛𝛼₂𝑥_𝑛⁻¹𝑥₂^𝑛−1∏_𝑖=1𝑥^𝛼𝑖_𝑖⋯𝑐₂𝛼_𝑛𝛼_𝑛−1𝑥_𝑛^𝑛−2∏_𝑖=1𝑥^𝛼𝑖_𝑖. (12)
Moreover, the determinant of ∇²𝑓₂(𝑥) can be computed and be shown by induction as det∇²𝑓₂(𝑥)=−𝑐₂^𝑛𝑛_𝑖=1𝛼_𝑖𝑥^𝑛𝛼𝑖_𝑖⁻²1−^𝑛_𝑖=1𝛼_𝑖. (13) We will complete the proof by discussing the following two cases.
Case i. If 1−∑^𝑛_𝑖=1𝛼_𝑖=0, then it is not hard to verify that ∇²𝑓₂(𝑥)𝑥 is a zero vector for any 𝑥∈ℝ^𝑛. Hence, 𝑥^𝑇∇²𝑓₂(𝑥)𝑥=0 which says that ∇²𝑓₂(𝑥) is a positive semidefinite matrix by definition. Therefore, 𝑓₂(𝑥) is convex under this case.
Case ii. If 1−∑^𝑛_𝑖=1𝛼_𝑖>0, then we know from (13) that Δ_𝑖=−𝑐₂^𝑖𝑖_𝑗=1𝛼_𝑗𝑥^𝑖𝛼𝑗_𝑗⁻²1−^𝑖_𝑗=1𝛼_𝑗, (14) where Δ_𝑖 denotes the 𝑖th principal minor of the Hessian matrix of 𝑓₂(𝑥). Note that 𝑐₂<0, 𝛼_𝑖>0 for all 𝑖=1,2,…,𝑛, and 1−∑^𝑛_𝑖=1𝛼_𝑖>0. Therefore, it can be seen that Δ_𝑖>0 for all 𝑖=1,2,…,𝑛, which implies (by Lemma 1(b)) that ∇²𝑓₂(𝑥) is a positive definite matrix. This says that 𝑓₂(𝑥) is strictly convex under this case.
From all the previous, the desired result follows.

We want to point out that our results also provide an alternative proof for the main result (Theorem7) of [28]. Indeed, Maranas and Floudas [28] further discuss another condition as follows: ∃𝑗suchthat𝛼_𝑗≥1−^𝑛_𝑖≠𝑗𝛼_𝑖,𝛼_𝑖≤0,∀𝑖≠𝑗,𝑖=1,2,…𝑛, (15) to guarantee that 𝑓₁ defined as in Proposition 3 is a convex function. Our approach can be also employed to verify this fact. To see this, we arrange all powers 𝛼_𝑖 in a decreasing order. In other words, without loss of generality, we assume that 𝛼₁>𝛼₂≥⋯≥𝛼_𝑛. (16) Notice that condition (15) implies that 𝛼₁ is positive and all the other 𝛼₂,…,𝛼_𝑛 are nonpositive with 𝛼₁≥1−∑^𝑛_𝑖=2𝛼_𝑖. As mentioned in Proposition 3, we only need to show that the function 𝑓₁(𝑥)=∏^𝑛_𝑖=1𝑥^𝛼𝑖_𝑖 is convex. By similar arguments as in the proof of Proposition 4, we know that Δ_𝑖=(−1)^𝑖𝑖_𝑗=1𝛼_𝑗𝑥^𝑖𝛼𝑗_𝑗⁻²1−^𝑖_𝑗=1𝛼_𝑗, (17) where Δ_𝑖 denotes the 𝑖th principal minor of the Hessian matrix of 𝑓₁(𝑥). From conditions (15) and (16), it is easily verified that (1−∑^𝑖_𝑗=1𝛼_𝑗)<0 for each 𝑖. It is also not hard to observe that ∏^𝑖_𝑗=1𝛼_𝑗 is positive if 𝑖 is odd and is negative if 𝑖 is even. In summary, there holds (−1)^𝑖𝑖_𝑗=1𝛼_𝑗𝑥^𝑖𝛼𝑗_𝑗⁻²<0. (18) The above two inequalities yield that Δ_𝑖>0 for each 𝑖. Thus, following the same arguments as in Proposition 4, we can conclude that 𝑓₁ is also a convex function under condition (15).

Remark 5. If 𝑐>0, 𝛼_𝑖<0, for 𝑖=1,2,…,𝑚, and 𝛼_𝑖>0, 𝑖=𝑚+1,𝑚+2,…,𝑛, then 𝑓(𝑥) can be converted into the following convex function: 𝑓^′(𝑥)=𝑐^𝑚_𝑖=1𝑥^𝛼𝑖_𝑖^𝑛_𝑖=𝑚+1𝑧^−𝛼𝑖_𝑖, (19) where 𝑧_𝑖=𝑥_𝑖⁻¹, 𝑖=𝑚+1,𝑚+2,…,𝑛.

Remark 6. If 𝑐<0, 𝛼_𝑖>0, for 𝑖=1,2,…,𝑚, and 𝛼_𝑖<0, 𝑖=𝑚+1,𝑚+2,…,𝑛, then 𝑓(𝑥) can be converted into the following convex function: 𝑓^′(𝑥)=𝑐^𝑚_𝑖=1𝑧^𝛼𝑖_𝑖^/𝑡𝑛_𝑖=𝑚+1𝑧^−𝛼𝑖_𝑖^/𝑡, (20) where (i)𝑡 is the smallest positive integer satisfying 1−(1/𝑡)(∑^𝑚_𝑖=1𝛼_𝑖−∑^𝑛_𝑖=𝑚+1𝛼_𝑖)≥0. (ii)𝑧_𝑖=𝑥^𝑡_𝑖, for 𝛼_𝑖>0, 𝑖=1,2,…,𝑚. (iii)𝑧_𝑖=𝑥_𝑖^−𝑡, for 𝛼_𝑖<0, 𝑖=𝑚+1,𝑚+2,…,𝑛.

Obviously, all signomial terms are converted into convex and concave terms by the proposed methods. All concave terms in the GP problem can be further represented by an effective piecewise linearization method as follows.

Conventional methods for linearizing a concave function with 𝑚 break points require 𝑚−1 binary variables. However, when 𝑚 becomes large, the computation will be very time consuming and may cause a heavy computational burden. An effective piecewise linearization method proposed by Huang [29] is presented in which only ⌈log₂(𝑚−1)⌉ binary variables are used.

Consider a grid of 𝑚 points throughout the interval of [𝑎₁,𝑎_𝑚] and denote the break points 𝑎₁,𝑎₂,…,𝑎_𝑚, where 𝑎₁<𝑎₂<⋯<𝑎_𝑚. 𝑓(𝑎_𝑘) is evaluated at each 𝑎_𝑘, and two adjacent break points (𝑎_𝑘,𝑓(𝑎_𝑘)) and (𝑎_𝑘+1,𝑓(𝑎_𝑘+1)) are connected with a straight line as shown in Figure 1.

Figure 1 

A piecewise linear function of 𝑓(𝑥).

Let 𝜃 be an integer, 1≤𝜃≤𝑚−1, expressed as 𝜃=1+^ℎ_𝑗=12^𝑗−1𝑢_𝑗, (21) where ℎ=⌈log₂(𝑚−1)⌉ and 𝑢_𝑗∈{0,1}. Given the values of 𝑢_𝑗, 1≤𝑗≤ℎ, there is a unique 𝜃.

Let 𝐺(𝜃)⊆{1,2,…,ℎ} be a set composed of all indices corresponding to a given value 𝜃, such that _{𝑗∈𝐺(𝜃)}2^𝑗−1=𝜃−1. (22)

For instance, 𝐺(1)=Ø, 𝐺(2)={1}, 𝐺(3)={2}, 𝐺(4)={1,2}, …, 𝐺(10)={1,4}, and so forth.

Denote by ‖𝐺(𝜃)‖ the number of elements in 𝐺(𝜃). For instance, given 1≤𝑥≤11, then ℎ=⌈log₂(𝑚−1)⌉=4. ‖𝐺(1)‖=0, ‖𝐺(2)‖=1, ‖𝐺(4)‖=2, and so forth. A list of all 𝜃, 𝑢_𝑗, 𝐺(𝜃), and ‖𝐺(𝜃)‖ is shown in Table 1.

Table 1 

Information of interval of 𝑥.

Now we intend to use 𝑘 as an index to indicate that the 𝑘th interval is activated, which results in 𝑎_𝑘≤𝑥≤𝑎_𝑘+1 and 𝑓(𝑥)=𝑓(𝑎_𝑘)+𝑠_𝑘(𝑥−𝑎_𝑘).

Define 𝑚−1 functions 𝑞_𝜃 based on the combination of 𝑢_𝑗 as follows: 𝑞_𝜃=‖𝐺(𝜃)‖−_{𝑢𝑗(𝜃)∈𝐷,𝑗∈𝐺}𝑢_𝑗+_{𝑢𝑗(𝜃)∈𝐷,𝑗∉𝐺}𝑢_𝑗, (23) where 𝜃=1,2,…,𝑚−1, 𝑗=1,2,…,ℎ, and ℎ=⌈log₂(𝑚−1)⌉.

The following proposition is deduced.

Proposition 7. Given a set of binary variables 𝐷={𝑢₁,𝑢₂,…,𝑢_ℎ}, where 𝑢_𝑗∈{0,1}, 𝑗=1,2,…,ℎ, and ℎ=⌈log₂(𝑚−1)⌉, there is one and only one 𝑘,1≤𝑘≤𝑚−1, such that (i)𝑞_𝑘=0 and ‖𝐺(𝑘)‖=∑_{𝑢𝑗∈𝐷,𝑗∈𝐺(𝑘)}𝑢_𝑗−∑_{𝑢𝑗∈𝐷,𝑗∉𝐺(𝑘)}𝑢_𝑗, (ii)𝑞_𝜃≥1 for 1≤𝜃≤𝑚−1, 𝜃≠𝑘.

Taking Table 1 for instance, given 1≤𝑥≤11, 𝑚=11. Let 𝐷={𝑢₁,𝑢₂,𝑢₃,𝑢₄}={1,0,0,1}; then we have 𝑞₁=‖𝐺(1)‖+𝑢₁+𝑢₂+𝑢₃+𝑢₄𝑞=0+1+0+0+1=2₂=‖𝐺(2)‖−𝑢₁+𝑢₂+𝑢₃+𝑢₄⋮𝑞=1−1+0+0+1=1₈=‖𝐺(8)‖−𝑢₁−𝑢₂−𝑢₃+𝑢₄𝑞=3−1−0−0+1=3₉=‖𝐺(9)‖+𝑢₁+𝑢₂+𝑢₃−𝑢₄𝑞=1+1+0+0−1=1₁₀=‖𝐺(10)‖−𝑢₁+𝑢₂+𝑢₃−𝑢₄=2−1+0+0−1=0. (24)

Since 𝑞₁₀=0 and 𝑞_𝜃≥1 for 1≤𝜃≤10, 𝜃≠10, then only 𝑎₁₀≤𝑥≤𝑎₁₁ is activated and results in 𝑓(𝑥)=𝑓(𝑎₁₀)+𝑠₁₀(𝑥−𝑎₁₀).

We then have the following theorem.

Theorem 8. Consider a piecewise linear function 𝑓(𝑥) of a concave function 𝑓(𝑥), where 𝑥 is divided into 𝑚−1 intervals by 𝑚 break points, 𝑎₁<𝑎₂<⋯<𝑎_𝑚. Given a set of binary variables 𝐷={𝑢₁,𝑢₂,…,𝑢_ℎ}, where ℎ=⌈log₂(𝑚−1)⌉,𝑥 and 𝑓(𝑥) can be expressed as the following linear inequalities: 𝑎_𝜃−𝑀𝑞_𝜃≤𝑥≤𝑎_𝜃+1+𝑀𝑞_𝜃,𝑓𝑎_𝜃+𝑠_𝜃𝑥−𝑎_𝜃−𝑀^′𝑞_𝜃≤𝑓(𝑥)≤𝑓𝑎_𝜃+𝑠_𝜃𝑥−𝑎_𝜃+𝑀^′𝑞_𝜃,𝑞_𝜃=‖𝐺(𝜃)‖−_{𝑢𝑗(𝜃)∈𝐷,𝑗∈𝐺}𝑢_𝑗+_{𝑢𝑗(𝜃)∈𝐷,𝑗∉𝐺}𝑢_𝑗, (25)(26)(27) where 𝑠_𝜃=(𝑓(𝑎_𝜃+1)−𝑓(𝑎_𝜃))/(𝑎_𝜃+1−𝑎_𝜃),𝜃=1,2,…,𝑚−1. 𝑀, and 𝑀^′ are large positive values,𝑀=max{𝑎_𝑚,1} and 𝑀^′=max{𝑓(𝑎₁),𝑓(𝑎₂),…,𝑓(𝑎_𝑚),1}.

Proof. Given the value of a specific 𝐷=𝐷⁰, 𝐷⁰={𝑢⁰₁,𝑢⁰₂,…,𝑢⁰_ℎ}, consider the following two cases.
Case i. There is one and only one integer 𝜃⁰, 𝜃⁰=1+∑^ℎ_𝑗=12^𝑗−1𝑢⁰_𝑗 such that (a)𝐺(𝜃⁰)={𝑗∣𝑢⁰_𝑗=1, (b)‖𝐺(𝜃⁰)‖=∑_{𝑗∈𝐺(𝜃}⁰₎𝑢⁰_𝑗, (c)∑_{𝑢0𝑗∈𝐷,𝑗∉𝐺(𝜃}⁰₎𝑢⁰_𝑗=0. Expressions (a), (b), and (c) force 𝑞_𝜃⁰=0 (from (27)). Then, 𝑎_𝜃⁰≤𝑥≤𝑎_𝜃⁰₊₁ is activated (from (25)) and results in 𝑓(𝑥)=𝑓(𝑎_𝜃⁰)+𝑠_𝜃⁰(𝑥−𝑎_𝜃⁰) (from (26)).
Case ii. For other 𝑚−1 integers 1≤𝜃^′≤𝑚, 𝜃^′≠𝜃⁰, since ‖𝐺(𝜃^′)‖>∑_{𝑗∈𝐺(𝜃′)}𝑢⁰_𝑗+∑_{𝑢0𝑗∈𝐷,𝑗∉𝐺(𝜃′)}𝑢⁰_𝑗, we then have 𝑞_𝜃′≥1. In this case, inequalities (25) and (26) are still correct but not being activated.

Compared with the conventional piecewise linearization methods [30–36], the number of newly added binary variables in the proposed method for a piecewise linear function with 𝑚 break points is significantly reduced from 𝑚−1 to ⌈log₂(𝑚−1)⌉. The algorithm is listed in Algorithm 1 based on Theorem 8.

Algorithm 1 (the proposed piecewise linearization method). Step 1. Identify the number of subintervals and choose the break points, 𝑎_𝑖, 𝑖=1,2,…,𝑚.
Step 2. Calculate the values of  ℎ=⌈log₂𝑚−1⌉ based on the number of break points selected and obtain 𝑞_𝜃, 𝜃=1,2,…,𝑚−1, based on the combination of 𝑢_ℎ as (27).
Step 3. Generate the set of constraints for 𝑥 and 𝑓(𝑥) as (25) and (26), respectively, and solve the resulting problem.

3. Numerical Examples

In this section, we have conducted some engineering design problems to evaluate the usefulness of the proposed methods.

Example 1. Consider the following engineering design problem of a speed reducer which was proposed by Golinski [37] as in Figure 2.

Figure 2 

Schematic of the welded beam problem, Example 2.

The objective of this problem is to minimize the weight of the speed reducer while satisfying a number of constraints imposed by gear and shaft design practices. There are seven design variables, 𝑥₁ (width of the gear face, cm), 𝑥₂ (teeth module, cm), 𝑥₃ (number of pinion teeth), 𝑥₄ (shaft 1 length between bearings, cm), 𝑥₅ (shaft 2 length between bearings, cm), 𝑥₆ (diameter of shaft 1, cm), and 𝑥₇ (diameter of shaft 2, cm). The constraints are characterized by ℎ₁ (upper bound on the bending stress of the gear tooth), ℎ₂ (upper bound on the contact stress of the gear tooth), ℎ₃, ℎ₄ (upper bounds on the transverse deflection of shafts 1, 2), ℎ₅, ℎ₆ (upper bounds on the stresses in shafts 1, 2), ℎ₇, ℎ₈, ℎ₉ (dimensional restrictions based on space and experience), and ℎ₁₀, ℎ₁₁ (dimensional requirements for shafts based on experience).

The problem can be formulated as min2.6179𝑥₁𝑥²₂𝑥²₃+11.7287𝑥₁𝑥²₂𝑥₃−33.84556𝑥₁𝑥²₂−1.508𝑥₁𝑥²₆−1.508𝑥₁𝑥²₇+7.4777𝑥³₆+7.4777𝑥³₇+0.7854𝑥₄𝑥²₆+0.7854𝑥₅𝑥²₇s.t.ℎ₁=−𝑥₁𝑥²₂𝑥₃ℎ+27.0≤0,₂=−𝑥₁𝑥²₂𝑥²₃ℎ+397.5≤0,₃=−𝑥₂𝑥₃𝑥⁴₆+1.93𝑥³₄ℎ≤0,₄=−𝑥₂𝑥₃𝑥⁴₇+1.93𝑥³₅ℎ≤0,₅=555025.0𝑥²₄+1.690×10⁷𝑥²₂𝑥²₃−12100.0𝑥²₂𝑥²₃𝑥⁶₆ℎ≤0,₆=555025.0𝑥²₅+1.575×10⁸𝑥²₂𝑥²₃−7225.0𝑥²₂𝑥²₃𝑥⁶₇ℎ≤0,₇=𝑥₂𝑥₃ℎ−40.0≤0,₈=−𝑥₁+5𝑥₂ℎ≤0,₉=𝑥₁−12𝑥₂ℎ≤0,₁₀=−𝑥₄+1.5𝑥₆ℎ+1.9≤0,₁₁=−𝑥₅+1.1𝑥₇+1.9≤0,2.6≤𝑥₁≤3.6,0.7≤𝑥₂≤0.8,17.0≤𝑥₃≤28.0,7.3≤𝑥₄≤8.3,7.3≤𝑥₅≤8.3,2.9≤𝑥₆≤3.9,5.0≤𝑥₇≤5.5. (28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40)(41)

This problem has been solved by Coello and Cortés [24] using genetic algorithm. The solution is (𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅,𝑥₆,𝑥₇) = (3.5, 0.7, 17.0, 7.300008, 7.715322, 3.350215, 5.286655) with the objective value 2994.3419.

For instance, the nonlinear terms 𝑥₁𝑥²₂𝑥²₃, 𝑥₁𝑥²₂𝑥₃, and 𝑥₄𝑥²₆ can be convexified as 𝑧⁻¹_1,1𝑧⁻²_2,1𝑧⁻²_3,1, 𝑧⁻¹_1,1𝑧⁻²_2,1𝑧⁻¹_3,1, and 𝑧⁻¹_4,1𝑧⁻²_6,1 by Remark 5, where 𝑧_𝑘,1=𝑥_𝑘⁻¹, 𝑘=1,2,…,7. The nonlinear terms −𝑥₁𝑥²₂, −𝑥₁𝑥²₆, and −𝑥₁𝑥²₇ can be convexified as −𝑧^1/3_1,2𝑧^2/3_2,2, −𝑧^1/3_1,2𝑧^2/3_6,2, and −𝑧^1/3_1,2𝑧^2/3_7,2 by Remark 6, where 𝑧_1,2=𝑥³₁, 𝑧_2,2=𝑥³₂, 𝑧_6,2=𝑥³₆, and 𝑧_7,2=𝑥³₇. The concave terms can be further represented by the proposed piecewise linearization method. Take 𝑧_1,1=𝑥₁⁻¹ for instance; suppose 41 equal-distance break points are used to linearize 𝑧_1,1=𝑥₁⁻¹, as given by 𝑎_1,1=2.600,𝑎_1,2=2.625,𝑎_1,3=2.650,…,and𝑎_1,41=3.600. Then 𝑧_1,1=𝑥₁⁻¹ can be represented by Algorithm 2 as follows.

Algorithm 2. Step 1. The proposed method requires the use of ⌈log₂(41−1)⌉=6 binary variables to linearize 𝑧_1,1=𝑥₁⁻¹ in a piecewise manner. Compared with the conventional piecewise linearization methods, the number of newly added binary variables for a piecewise linear function with 41 break points is significantly reduced from 40 to 6.
Step 2. ℎ_1,1=⌈log₂(41−1)⌉=6. A set of 𝑞_1,𝜃, 𝜃=1,2,…,40, based on the combination of 𝑢_1,ℎ1,1 is expressed as (42)𝑞_1,1=0+𝑢_1,1+𝑢_1,2+𝑢_1,3+𝑢_1,4+𝑢_1,5+𝑢_1,6𝑞_1,2=1−𝑢_1,1+𝑢_1,2+𝑢_1,3+𝑢_1,4+𝑢_1,5+𝑢_1,6⋮𝑞_1,40=4−𝑢_1,1−𝑢_1,2−𝑢_1,3+𝑢_1,4+𝑢_1,5−𝑢_1,6. (42)
Step 3. Based on Theorem 8, 𝑥₁ and 𝑧_1,1 can be expressed as the following linear inequalities:2.600−𝑀𝑞_1,1≤𝑥₁≤2.625+𝑀𝑞_1,12.625−𝑀𝑞_1,2≤𝑥₁≤2.650+𝑀𝑞_1,2⋮3.575−𝑀𝑞_1,40≤𝑥₁≤3.600+𝑀𝑞_1,40,2.600⁻¹+𝑠_1,1𝑥₁−2.600−𝑀^′𝑞_1,1≤𝑧_1,1≤2.600⁻¹+𝑠_1,1𝑥₁−2.600+𝑀^′𝑞_1,12.625⁻¹+𝑠_1,2𝑥₁−2.625−𝑀^′𝑞_1,2≤𝑧_1,1≤2.625⁻¹+𝑠_1,2𝑥₁−2.625+𝑀^′𝑞_1,2⋮3.575⁻¹+𝑠_1,40𝑥₁−3.575−𝑀^′𝑞_1,40≤𝑧_1,1≤3.575⁻¹+𝑠_1,40𝑥₁−3.575+𝑀^′𝑞_1,40, (43)(44) where 𝑠_1,𝜃=(𝑎⁻¹_1,𝜃+1−𝑎⁻¹_1,𝜃)/(𝑎_1,𝜃+1−𝑎_1,𝜃), 𝜃=1,2,…,40. 𝑀 and 𝑀^′ are large positive values.

Solving the transformed program with MATLAB R13, Table 2 shows the results computed by Coello and Cortés [24] and the proposed method.

Table 2 

Comparison of Example 1.

Example 2 (see [38]). Consider the design of a welded beam which is to be rigidly attached to a fixed support and is subjected to a load as illustrated in Figure 3. The welded beam is to consist of 1010 steel and is to support a force 𝐹 of 6000 lb. The length 𝐿 is assumed to be specified at 14 in.

Figure 3 

Golinski’s speed reducer problem, Example 1.

Three methods (brazed, bolted, and welded connections) of attachment and three materials (steel, brass, and aluminum) are considered. The choices of attachment and materials are to be made, which minimizes cost while satisfying a limit on the allowable bending stress, the buckling load, and the allowable deflection at the load application point. The various materials have different elastic, shear moduli, and different costs.

The design variables are defined as 𝑥₁ (a continuous variable for thickness of weld), 𝑥₂ (a continuous variable for width of beam), 𝑥₃ (a continuous variable for length of connection), 𝑥₄ (a continuous variable for thickness of beam), 𝑥₅ (a continuous variable for diameter of bolts), 𝑥₆ (a 0-1 variable for brazed connections), 𝑥₇ (a 0-1 variable for bolted connections), 𝑥₈ (a 0-1 variable for welded connections), 𝑥₉ (a 0-1 variable for steel material), 𝑥₁₀ (a 0-1 variable for brass material), and 𝑥₁₁ (a 0-1 variable for aluminum material).

The problem can be formulated as (45)–(50) min1.2𝑥₂𝑥₃𝑥₆+𝑥₅𝑥₇+𝑥²₁𝑥₃𝑥₈+0.04811𝐿𝑥₂𝑥₄𝑥₉+0.04811𝑥₂𝑥₃𝑥₄𝑥₉+0.06𝐿𝑥₂𝑥₄𝑥₁₀+0.06𝑥₂𝑥₃𝑥₄𝑥₁₀+0.015𝐿𝑥₂𝑥₄𝑥₁₁+0.015𝑥₂𝑥₃𝑥₄𝑥₁₁s.tℎ₁=−𝜎_allowable+6𝐹𝐿𝑥²₂𝑥₄ℎ≤0,₂=−𝑃_𝑐ℎ+6000≤0,₃=−𝛿_max+4𝐹𝐿³𝐸𝑥³₂𝑥₄𝑥≤0,₆+𝑥₇+𝑥₈𝑥=1,₉+𝑥₁₀+𝑥₁₁=1, (45)(46)(47)(48)(49)(50) where 𝜎_allowable is the maximum bar bending stress. 𝑃_𝑐 is the bar buckling load which is a function of the beam dimensions and material properties. 𝛿_max is the maximum bar deflection. 𝐸 is the Young’s modulus, 30×10⁶ psi. More detailed description including the derivation of 𝑃_𝑐 can be found in Reklaitis et al. [38].

Following the same procedure as in Example 1, the nonlinear terms 𝑥₂𝑥₃𝑥₆, 𝑥₅𝑥₇, 𝑥²₁𝑥₃𝑥₈, 𝑥₂𝑥₄𝑥₉, 𝑥₂𝑥₃𝑥₄𝑥₉, 𝑥₂𝑥₄𝑥₁₀, 𝑥₂𝑥₃𝑥₄𝑥₁₀, 𝑥₂𝑥₄𝑥₁₁, and 𝑥₂𝑥₃𝑥₄𝑥₁₁ can be convexified by Remark 5. The concave terms are further represented by the proposed piecewise linearization method. Solving the transformed program with MATLAB R13, we can obtain the global optimum with (𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅,𝑥₆,𝑥₇,𝑥₈,𝑥₉,𝑥₁₀,𝑥₁₁)=(0.320,14.109,5.720,0.321,4.992,0,0,1,0,0,1) and the objective value is 1.969 within the tolerance error 0.0001.

4. Conclusions

This paper proposes an optimization approach for solving geometric programming problems. Our approach is first to convert all signomial terms in GP into convex and concave terms by the proposed methods. Then the concave terms are further treated with the proposed piecewise linearization method where only ⌈log₂(𝑚−1)⌉ binary variables are used. It has the following features.(i)It offers more convenient and efficient means of expressing a piecewise linear function. (ii)Fewer 0-1 variables are used. (iii)The computational results show that the proposed method is much more efficient and faster than the conventional one, especially when the number of break points becomes large.

Numerical examples in real applications are illustrated to demonstrate the usefulness of the proposed methods.

Acknowledgments

The author would like to thank Distinguished Chair Professor Gwo-Hshiung Tzeng for reading the first draft of this paper and giving him many useful comments to improve this paper. The author also thanks the anonymous referees for their careful reading and fruitful comments on the paper. The author especially acknowledges the National Science Council, Taiwan, for financially supporting this research under Grant no. NSC 102–2410-H-141-012-MY2.