Convex restrictions in physical design.

In a physical design problem, the designer chooses values of some physical parameters, within limits, to optimize the resulting field. We focus on the specific case in which each physical design parameter is the ratio of two field variables. This form occurs for photonic design with real scalar fields, diffusion-type systems, and others. We show that such problems can be reduced to a convex optimization problem, and therefore efficiently solved globally, given the sign of an optimal field at every point. This observation suggests a heuristic, in which the signs of the field are iteratively updated. This heuristic appears to have good practical performance on diffusion-type problems (including thermal design and resistive circuit design) and some control problems, while exhibiting moderate performance on photonic design problems. We also show in many practical cases there exist globally optimal designs whose design parameters are maximized or minimized at each point in the domain, i.e., that there is a discrete globally optimal structure.

Introduction

Computer-aided physical design has become an important tool in many fields including photonics[1,2], mechanical design[3], circuit design[4,5], and thermal design[6,7]. In many cases, the design problem is formulated as a constrained nonconvex optimization problem which is then approximately minimized using local optimization methods such as ADMM[8], evolutionary algorithms[5], and the method of moving asymptotes[6], among many others.

More generally, a physical design problem can be phrased in the following way: we are allowed to choose some design parameters (e.g., the permittivity in photonic design or the conductances in diffusion design) at each point in a domain, within some limits, in order to minimize an objective function of the field (this can be, e.g., the electric field in photonic design, or a vector containing the potentials, flows, and potential differences in diffusion design). The constraints specify the physics of the problem, connecting the design variables to the field variables (e.g., Maxwell’s equations in photonics, or a diffusion equation such as the heat equation in diffusion design). We note that, in many cases, the physics constraints are linear equations in the field variables (when the design parameters are held constant), and linear equations in the design parameters (when the fields are held constant), which has led to some heuristics with good performance[8].

There has been recent interest in understanding global properties of solutions for physical design problems: lower bounds for optimal design objectives in photonic design have been studied via the use of convex relaxations found by physical arguments[9,10], duality theory[11-13], among others[14]. We instead analyze a convex restriction (see[15], Sect. 2.1) of the physical design problem, potentially providing another approach for analyzing properties of global solutions and for creating fast heuristics.

In this paper, we consider a simple (but very general) formulation of a class of physical design problems which includes problems in thermal design, photonic inverse design with scalar fields and convex objectives, and some types of control problems. This formulation offers some insights into the properties of global solutions for these problems. For example, in many practical cases, problems with linear objectives can be shown to have optimal extremal designs (in the case of physical design) or bang-bang controls (in the case of control). As another example, we observe that it suffices to know only the sign of a subset of variables in order to globally solve the problem efficiently, even though the original problem is NP-hard. The formulation also suggests a heuristic which appears to have good performance for many kinds of physical design problems, and we give numerical examples of this heuristic applied to a few different problems.

General problem formulation

We consider a problem of the formwhere is a convex function over our variables and , is a convex constraint set, and is our design variable whose limits are . While apparently simple, many physical design problems can be expressed as instances of problem (1); we show a few examples in the “Applications” section. We call (x, u, v) the field (corresponding to, e.g., the electric field in photonic design) and the design parameters (corresponding to, e.g., the permittivity in photonic design). We say that is extremal whenever for each . The physics of the problem is encoded in the constraints and .

In this problem, the convex set can be any convex set specifying constraints on the variables (x, u, v), such as linear equality constraints. On the other hand, the design parameters enter in a very specific way: as a diagonal term relating u and v. Another way to say this is that each design parameter is the ratio of two field parameters, and .

We note that the problem (1) is convex in (x, u, v) whenever is fixed, and convex in whenever v is fixed. In practice, there has been great success in applying heuristics for approximately minimizing instances of (1) using this observation[16].

Absolute upper bound formulation

Problem (1) is equivalent towhere the absolute value is taken elementwise. The variables of problem (2) are and , while and are constants. Note that is the middle value of the physical parameter interval, and is the radius, i.e., half the range or width of the interval.

The equivalence between problems (1) and (2) can be seen by noting that, for every feasible (x, u, v, w) for problem (2) we can set,for . Then, is feasible for (1), with the same objective value. Note that, if , any choice of would suffice.

Similarly, for any that is feasible for (1), we can setand then (x, u, v, w) is feasible for problem (2) with the same objective value.

We will refer to problem (2) as the absolute-upper-bound formulation of problem (1). This problem, like problem (1), is nonconvex due to the inequality , and is hard to solve exactly.

NP-hardness

We can reduce any mixed-integer convex program (MICP) to an instance of (2), implying that this problem is hard, as any instance of an NP-complete problem is easily reducible to instances of the MICP problem[17].

The reduction follows since we can force v to be binary in problem (2). First, choose , (and therefore ), and add to the constraint set. This immediately implies that . Adding the convex constraint to the constraint set , yields , as required. Since and f can be otherwise freely chosen, the result follows.

Known signs

If the signs of an optimal are known for problem (2), then the problem becomes convex. We can see this as follows. If is known, then we can solve the following convex problem[18], Sect. 4:where is the elementwise product of s and v. Note that (and its associated values of , , and ) are feasible for this instance of (4) since , which implies that a solution of this instance of (4) must be globally optimal for (2).

Global solution

Note that problem (4) generates a family of optimization problems over the set of possible signs, . This suggests a simple, if inefficient, way to globally solve problem (2) and therefore problem (1): solve problem (4) for the possible signs, , to obtain optimal values for each set of signs s. A solution for any optimal set of signs, , is then a solution to (2) and therefore to (1).

Of course, this algorithm may not be useful in practice for anything but the smallest values of m, but it implies that solving problem (1) requires solving only a finite number of convex problems.

Extremality principle

The rewriting given in (4) also yields an interesting insight. If problem (4) is a feasible linear program and is an affine set with , i.e., for each there exists a and an such that , then there exists a solution of (4) such that all entries of the inequality hold at equality (see, e.g.,[19], Sect. 2.6). This rewriting then implies that there exists an optimal design for which is extremal, by (3). A numerical example of this principle is found in the “Thermal design” section.

Sign flip descent

Since problem (4) generates a family of optimization problems parametrized by the sign vector , we can view the original physical design problem (1) as a problem of choosing an optimal Boolean vector. A simple way of approximately optimizing (2) is: at each iteration i, start with some sign vector and solve (4) to obtain an optimal value . We then consider a rule for proposing a new sign vector, say , for which we again solve (4) and then obtain a new optimal value . If , we then keep this new sign vector, i.e., we set , and repeat the procedure; otherwise, we discard by setting , and repeat the procedure, proposing a new sign vector in the next iteration. This is outlined in algorithm 1.

By construction, any algorithm of the form of algorithm 1 is a descent algorithm since each iteration is feasible and the objective value is decreasing on each iteration. We outline two possible rules for proposing new sets of signs at each iteration.

Greedy sign rule

A simple rule for choosing signs is to begin at iteration k with some set of signs . We then define a new set of signs with except at the kth entry where we have (or, if then we pick the entry at index , i.e., such that the entries are changed, one-by-one, in a round-robin fashion). We stop whenever flipping any one entry of does not yield a lower objective value.

The greedy sign rule has two useful properties. First, the rule guarantees local optimality in the following sense: if algorithm 1 returns , then changing any one sign of will not decrease the objective value. Second, the rule terminates in finite time, since the corresponding algorithm is a descent algorithm and there are a finite number of possible sign vectors. On the other hand, the algorithm is often slow for anything but the smallest designs: to reach the terminating condition, we have to solve at least m convex optimization problems.

Field-based rule

Another simple rule that appears to work very well in practice is based on the observation that, for many choices of sign vectors , the inequality has many entries of v that are zero. If is zero for some index , this suggests that the sign might have been originally set incorrectly: in this case, we propose a new vector which is equal at all entries for which is nonzero and has opposite sign at all entries i for which is zero.

Note that this new proposed vector will always have an optimal value which is at least as small as the optimal value for , i.e., . This observation, coupled with the proposed rule, suggests that we should stop whenever there are no signs left to flip, or whenever the iterations stop decreasing as quickly as desired, i.e., whenever .

While this rule does not necessarily guarantee local optimality, it always terminates in finite time with the given stopping conditions and appears to work well in practice (requiring, in comparison to the greedy sign rule, much fewer than m iterations to terminate) as shown in the “Numerical examples” section.

Applications

We describe a few interesting design problems that reduce to problems of the form of (1).

Diagonal physical design

As in, e.g.,[11], many physical design problems can be written in the following way:where describes the physics of the problem, while describes the excitation, and are the design parameters of the system, chosen to minimize some convex objective function of the field . Our variables in this problem are the field z and the design parameters .

We can write a problem of the form of (5) as a problem of the form (1) by introducing a new variable u with constraint and rewriting the equality constraint of (5) with this new variable, . As the set of (z, u) satisfying forms a convex (in fact, affine) set, the resulting problem,is of the form of (1) which can be easily rewritten into the form of (2).

Static diffusion design

Consider a flow problem on a graph where we choose the conductance across each edge , constrained to satisfy , to minimize some function of the potentials , given some sources .

To compactly write the conditions this system must satisfy, let the matrix be the incidence matrix for the graph G defined to be (see[20], Sect. 7.3):

We can then write the steady-state diffusion equation aswhere can be recognized as the graph Laplacian of G with edge weights g. This equation can also be seen as the discrete form of the heat equation on a graph G[21].

The corresponding optimization problem is then an instance of (1):where we have introduced two new variables , in addition to the potential and the conductances . As before, is the incidence matrix, are the sources at each node, while is a vector such that is the average temperature over the desired region.

Dynamic diffusion control

Similarly to the “Static diffusion design” section, we can consider the time-varying generalization of (6) given byat each time , with step size . Here, is the heat capacity of each node and , while are the inputs given to the system, is a matrix mapping these inputs to the power added or removed from each node, are the conductances at each node, and is the temperature at each node.

In this case, we can minimize any convex function of the temperatures and inputs by appropriately choosing the conductances and inputs:where, as before, we have introduced the variables , for each and .

We can see problem (8) as a nontraditional control problem. A particular example is: we have a set of rooms with temperatures at time t which we wish to keep within some comfortable temperature range. We are allowed to open and close vents (equivalently, change the conductances at each time t) and turn on and off heat pumps (via the control variable ), while paying a cost for the latter. A simple question could be: what is an optimal set of actions such that the input cost is minimized while keeping the temperatures within some specified bounds? We show a simple example of this in the “Temperature control” section.

Numerical examples

Julia[22] code for all examples in this section is available in the following Github repository: angeris/pd-heuristic. We use the JuMP modeling language[23] to interface with Mosek[24]. All times reported are on a 2015 2.9 GHz dual-core MacBook Pro.

Photonic design

In this example, we wish to choose the speed of a wave satisfying Helmholtz’s equation at each point in some domain in order to minimize a convex function of the field.

Helmholtz’s equation

More specifically, the speed of the wave is chosen such that the field at a specific frequency with excitation satisfies Helmholtz’s equation,at each point . Additionally, we require that the chosen speeds are bounded such that at each point , and we assume Dirichlet boundary conditions such that for , i.e., we require the field to be zero at every point on the boundary of the domain. In electromagnetics, this condition corresponds to having a perfect conductor at the boundary.

In this case (as in[11], Sect. 5.1), we will work with a discretized form of (9) where is the discretized field (), is the discretized excitation (), is the discretized speed of the wave (c), and is the discretized version of the Laplacian operator (), such thatapproximates (9) at each point for . We assume that the discretization is such that is a box.

Problem data

In this case, the problem data are given by , with , while the convex objective function is given bywhere is the box indicated in Fig. 1, and the excitation b is defined asfor each , where is the box indicated in Fig. 1. Here, and . We set the tolerance parameter of the algorithm to . We initialize the algorithm by finding a solution to Eq. (10) with and use the signs of this solution as the initial sign vector.

Figure 1

Approximately optimal photonic design. The leftmost figure specifies in purple, the center specifies , while the rightmost figure gives the design, .

Numerical results

With the given problem data, the algorithm terminates at 102 iterations with a total time of about 4 minutes, roughly around 2 seconds per iteration. This time could be very much shortened since the current implementation does not warm-start any of the current iterations, essentially solving the problem from scratch at each iteration. The final design is shown in Fig. 1 and its final field is shown in Fig. 2.

Figure 2

Field for approximately optimal design.

Thermal design

In this design problem, as in the “Static diffusion design” section, we seek to set the conductances on a graph in order to minimize the average temperature of a subset of points in the center of a 2D grid of size , given a heat source and a heat sink at opposite corners of the 2D grid. This is an instance of the diffusion problem where is the incidence matrix of the grid and are the heat sources and sinks. This problem can be written as an instance of (7) where the potentials are the temperatures at each point in the grid.

Our convex objective function is given bywhere is a vector such that if vertex i lies in the center square of size while otherwise. There is a heat source set at the bottom left corner of the grid and a heat sink set at the top right corner of the grid. We set the minimal and maximal conductances as and at each edge.

We approximately optimize the conductances in this problem by using the field-based heuristic described in the “Sign flip descent” section. The directions are initialized by solving the problem with uniform conductances.

A small example is given in Fig. 3 with (which shows the chosen directions of flow), while a relatively large design is given in Fig. 4 with . In both figures, thick edges indicate that conductance is maximized at that edge while thin edges indicate that conductance is minimized (see the extremality principle in the “General problem formulation” section for more details). The color of each node indicates the potential value, with red values indicating a higher potential and blue values indicating a lower one. We note that our heuristic recovers similar tendril-like patterns to those found in, e.g.[7], Sect. 4 .

Figure 3

Approximately optimal design for . Arrows indicate the direction of flow used for this design, colors indicate the temperature at each node, while edge thickness indicates the conductance at each edge. The grey box indicates the center square.

Figure 4

Approximately optimal design for .

With the provided data, the heuristic terminates after 7 iterations, taking a total time of around .4 seconds in the case with , with an objective value of about .115. The case with terminates after 14 iterations, taking a total time of around 20.5 seconds with an objective value of approximately .239.

Temperature control

In this example, we wish to keep the temperature of two rooms in a range of desired temperatures by appropriately closing and opening vents to the outside and between rooms and turning heat pumps on and off at specified times, while minimizing the total power consumption. We will also require that the controls and the temperatures be periodic.

We can write this as an instance of problem (8) withand A is the incidence matrix of the graph shown in Fig. 5, whilewhere . Since we will require that the room temperatures be periodic, we then have

Figure 5

Graph set up for the temperature control problem. Here, is the ambient temperature at time t, while and are the temperatures of rooms 1 and 2, respectively. The are the conductances of the indicated edges.

Finally, we will require that the temperatures remain in some a range,while minimizingwhere and is a small regularization parameter that ensures the resulting trajectories are smooth.

We initialize the problem with the signs given by assuming that for all and using the heat pumps to ensure the temperature in the rooms remains above 65 and below 75.

We approximately optimize this instance using the field-based heuristic outlined in “Sign flip descent”, with the result shown in Fig. 6. With the provided data, the heuristic terminates in 3 iterations, with a total time of around 1.56 s. The final approximately optimized problem has an objective value of around 836.

Figure 6

Approximately optimal control.

Conclusion

This paper presented a new problem formulation and an associated heuristic which may be of practical use for a general class of physical design problems, which appears to have good practical performance on many different kinds of physical design problems. Additionally, this problem formulation implies a few interesting facts, most notably that the class of problems can be efficiently solved even when only the signs of an optimal solution are known and that, in a few important cases, there exist globally optimal extremal designs.

Future work

There are several notable exceptions to the class of problems which are included in the formulation given in (1), with the most important being designs whose parameters are constrained to be equal. This means that, at the moment, a direct application to photonic design in three dimensions, the usual photonic design problem with complex fields, circuit design with complex impedances, or multi-scenario physical design, is not possible with the current problem formulation. We suspect a suitable generalization of (1) might yield similarly interesting insights and, potentially, new heuristics for physical design.