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.

www.nature.com/scientificreports/ where f : R n × R m × R m → R is a convex function over our variables x ∈ R m and u, v ∈ R n , C ⊆ R n × R m × R m is a convex constraint set, and θ ∈ R n is our design variable whose limits are θ min , θ max ∈ R n . 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 θ i ∈ {θ min i , θ max i } for each i = 1, . . . , m . The physics of the problem is encoded in the constraints (x, u, v) ∈ C and u = diag(θ)v.
In this problem, the convex set C 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 θ i is the ratio of two field parameters, u i and v i .
We note that the problem (1) is convex in (x, u, v) whenever θ is fixed, and convex in (x, u, θ) 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 to where the absolute value is taken elementwise. The variables of problem (2) are x ∈ R m and u, v, w ∈ R n , while θ = (θ max + θ min )/2 and ρ = (θ max − θ min )/2 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 i = 1, . . . , m . Then, (x, u, v, θ) is feasible for (1), with the same objective value. Note that, if v i = 0 , any choice Similarly, for any (x, u, v, θ) that is feasible for (1), we can set and 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 |w| ≤ |v| , 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 θ = 0 , ρ = 1 (and therefore u = w ), and add u = 1 to the constraint set. This immediately implies that 1 ≤ |v| . Adding the convex constraint |v| ≤ 1 to the constraint set C , yields v ∈ {±1} n , as required. Since C and f can be otherwise freely chosen, the result follows.
Known signs. If the signs of an optimal v ⋆ are known for problem (2), then the problem becomes convex.
We can see this as follows. If s = sign(v ⋆ ) ∈ {±1} m is known, then we can solve the following convex problem 18 ,Sect. 4: where s • v is the elementwise product of s and v. Note that v ⋆ (and its associated values of x ⋆ , u ⋆ , and w ⋆ ) are feasible for this instance of (4) since |v ⋆ | = s • v ⋆ , 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, s ∈ {±1} m . This suggests a simple, if inefficient, way to globally solve problem (2) and therefore www.nature.com/scientificreports/ problem (1): solve problem (4) for the 2 m possible signs, s ∈ {±1} m , to obtain optimal values p ⋆ (s) for each set of signs s. A solution (x ⋆ , u ⋆ , v ⋆ , w ⋆ ) for any optimal set of signs, s ⋆ ∈ argmin s∈{±1} m p ⋆ (s) , 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 C is an affine set with {u | (x, u, v) ∈ C } = R m , i.e., for each u ∈ R m there exists a v ∈ R m and an x ∈ R n such that (x, u, v) ∈ C , then there exists a solution of (4) such that all entries of the inequality |w| ≤ s • v 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 s ∈ {±1} m , 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 s i ∈ {±1} m and solve (4) to obtain an optimal value p i . We then consider a rule for proposing a new sign vector, say s i ∈ {±1} m , for which we again solve (4) and then obtain a new optimal value p i . If p i < p i , we then keep this new sign vector, i.e., we set s i+1 =s i , and repeat the procedure; otherwise, we discard s i by setting s i+1 = s i , 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 s k . We then define a new set of signs s k with s k = s k except at the kth entry where we have s k k = −s k k (or, if k > m then we pick the entry at index 1 + (k − 1 mod m) , i.e., such that the entries are changed, one-by-one, in a round-robin fashion). We stop whenever flipping any one entry of s k 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 s ⋆ , then changing any one sign of s ⋆ 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 s k , the inequality |w| ≤ s k • v has many entries of v that are zero. If v i is zero for some index i = 1, . . . , m , this suggests that the sign s k i might have been originally set incorrectly: in this case, we propose a new vector s k which is equal s k at all entries i = 1, . . . , m for which v i is nonzero and has opposite sign at all entries i for which v i is zero.
Note that this new proposed vector will always have an optimal value p k which is at least as small as the optimal value for s k , i.e., p k ≤ p k . 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 p k − p k+1 < ε.
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 A ∈ R n×n describes the physics of the problem, while b ∈ R n describes the excitation, and θ ∈ R n are the design parameters of the system, chosen to minimize some convex objective function f : R n → R of the field z ∈ R n . 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 u = diag(θ)z and rewriting the equality constraint of (5) with this new variable, Az + u = b . As the set of (z, u) satisfying Az + u = b 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 G = (V , E) where we choose the conductance g k ∈ R across each edge k ∈ E , constrained to satisfy g min k ≤ g k ≤ g max k , to minimize some function f : R |V | → R of the potentials e ∈ R |V | , given some sources s ∈ R |V | .
To compactly write the conditions this system must satisfy, let the matrix A ∈ R |V |×|E| 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 as where Adiag(g)A T 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 w, v ∈ R |E| , in addition to the potential e ∈ R |V | and the conductances g ∈ R |E| . As before, A ∈ R |V |×|E| is the incidence matrix, s ∈ R |V | are the sources at each node, while c ∈ R n is a vector such that c T e is the average temperature over the desired region.

Dynamic diffusion control.
Similarly to the "Static diffusion design" section, we can consider the timevarying generalization of (6) given by at each time t = 1, . . . , T , with step size h > 0 . Here, c ∈ R |V | ++ is the heat capacity of each node and C = diag(c) , while u t ∈ R n are the inputs given to the system, B ∈ R |V |×n is a matrix mapping these inputs to the power added or removed from each node, g t ∈ R |V | are the conductances at each node, and e t ∈ R |V | 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 v t , w t ∈ R |E| , for each t ∈ [T] and [T] = {1, . . . , T}.
A ij = +1 edge j points to node i −1 edge j points from node i 0 otherwise.
minimize f (e) subject to v = A T e Aw = s w = diag(g)v g min ≤ g ≤ g max ,  (8) as a nontraditional control problem. A particular example is: we have a set of rooms with temperatures e t 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 g t at each time t) and turn on and off heat pumps (via the control variable u t ), 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 e t within some specified bounds? We show a simple example of this in the "Temperature control" section. Helmholtz's equation. More specifically, the speed of the wave c : → R ++ is chosen such that the field ψ : � → R at a specific frequency ω ∈ R + with excitation φ : � → R satisfies Helmholtz's equation, at each point (x, y) ∈ � . Additionally, we require that the chosen speeds are bounded such that 0 < c min (x, y) ≤ c(x, y) ≤ c max (x, y) at each point (x, y) ∈ � , and we assume Dirichlet boundary conditions such that ψ(x, y) = 0 for (x, y) ∈ ∂� , 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.

Numerical examples
In this case (as in 11 , Sect. 5.1), we will work with a discretized form of (9) where z ∈ R n is the discretized field ( ψ ), b ∈ R n is the discretized excitation ( φ ), θ ∈ R n is the discretized speed of the wave (c), and A ∈ R n×n is the discretized version of the Laplacian operator ( ∇ 2 ), such that approximates (9) at each point (x i , y i ) ∈ � for i = 1, . . . , n . We assume that the discretization is such that is a 1 × 1 box.
Problem data. In this case, the problem data are given by ω = 4π , with n = 101 × 101 = 10201 , while the convex objective function f : R n → R is given by where B ⊆ {1, . . . , n} is the box indicated in Fig. 1, and the excitation b is defined as for each i = 1, . . . , n , where S ⊆ {1, . . . , n} is the box indicated in Fig. 1. Here, θ min = 1 and θ max = 2 . We set the tolerance parameter of the algorithm to ε = 10 −4 . We initialize the algorithm by finding a solution to Eq. (10) with θ = (θ max + θ min )/2 and use the signs of this solution as the initial sign vector. www.nature.com/scientificreports/ 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.
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 m × m , given a heat source and a heat sink at opposite corners of the 2D grid. This is an instance of the diffusion problem where A ∈ R |E|×|V | is the incidence matrix of the grid and s ∈ R |V | are the heat sources and sinks. This problem can be written as an instance of (7) where the potentials e ∈ R |V | are the temperatures at each point in the grid.
Problem data. Our convex objective function f : R |V | → R is given by where c ∈ R |V | is a vector such that c i = 1 if vertex i lies in the center square of size ⌊(m − 1)/4⌋ × ⌊(m − 1)/4⌋ while c i = 0 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 g min = 1 and g max = 10 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.
Numerical results. A small example is given in Fig. 3 with m = 11 (which shows the chosen directions of flow), while a relatively large design is given in Fig. 4 with m = 51 . 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  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.
Problem data. We can write this as an instance of problem (8) with and A is the incidence matrix of the graph shown in    www.nature.com/scientificreports/ where T = 300 . Since we will require that the room temperatures be periodic, we then have Finally, we will require that the temperatures remain in some a range,   where h = 1/T and η = 10 −4 is a small regularization parameter that ensures the resulting trajectories are smooth.
We initialize the problem with the signs given by assuming that g t = (g min + g max )/2 for all t = 1, . . . , T − 1 and using the heat pumps u t to ensure the temperature in the rooms remains above 65 and below 75.
Numerical results. 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.

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.