Original by Lieven Vanderberghe.
Adapted to CVX by Argyris Zymnis, 12/4/2005.
Modified by Michael Grant, 3/8/2006.
Adapted to CVXPY, with cosmetic modifications, by Judson Wilson, 5/26/2014.
Topic References:
We consider the problem of sizing a clock mesh, so as to minimize the total dissipated power under a constraint on the dominant time constant. The numbers of nodes in the mesh is $N$ per row or column (thus $n=(N+1)^2$ in total). We divide the wire into m segments of width $x_i$, $i = 1,\dots,m$ which is constrained as $0 \le x_i \le W_{\mbox{max}}$. We use a pi-model of each wire segment, with capacitance $\beta_i x_i$ and conductance $\alpha_i x_i$. Defining $C(x) = C_0+x_1 C_1 + x_2 C_ 2 + \cdots + x_m C_m$ we have that the dissipated power is equal to $\mathbf{1}^T C(x) \mathbf{1}$. Thus to minimize the dissipated power subject to a constraint in the widths and a constraint in the dominant time constant, we solve the SDP \begin{array}{ll} \mbox{minimize} & \mathbf{1}^T C(x) \mathbf{1} \ \mbox{subject to} & T_{\mbox{max}} G(x) - C(x) \succeq 0 \ & 0 \le xi \le W{\mbox{max}}. \end{array}
import cvxpy as cvx
import numpy as np
import scipy as scipy
import matplotlib.pyplot as plt
# Show plots inline in ipython.
%matplotlib inline
# Plot properties.
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
font = {'family' : 'normal',
'weight' : 'normal',
'size' : 16}
plt.rc('font', **font)
# Computes the step response of a linear system.
def simple_step(A, B, DT, N):
n = A.shape[0]
Ad = scipy.linalg.expm((A * DT))
Bd = (Ad - np.eye(n)) * B
Bd = np.linalg.solve(A, Bd)
X = np.mat(np.zeros((n, N)))
for k in range(1, N):
X[:, k] = Ad * X[:, k-1] + Bd;
return X
#
# Circuit parameters.
#
dim=4 # Grid is dimxdim (assume dim is even).
n=(dim+1)**2 # Number of nodes.
m=2*dim*(dim+1) # Number of wires.
# 0 ... dim(dim+1)-1 are horizontal segments
# (numbered rowwise);
# dim(dim+1) ... 2*dim(dim+1)-1 are vertical
# (numbered columnwise)
beta = 0.5 # Capacitance per segment is twice beta times xi.
alpha = 1 # Conductance per segment is alpha times xi.
G0 = 1 # Source conductance.
C0 = np.array([ ( 10, 2, 7, 5, 3),
( 8, 3, 9, 5, 5),
( 1, 8, 4, 9, 3),
( 7, 3, 6, 8, 2),
( 5, 2, 1, 9, 10) ])
wmax = 1 # Upper bound on x.
#
# Build capacitance and conductance matrices.
#
CC = np.zeros((dim+1, dim+1, dim+1, dim+1, m+1))
GG = np.zeros((dim+1, dim+1, dim+1, dim+1, m+1))
# Constant terms.
# - Reshape order='F' is fortran order to match original
# version in MATLAB code.
CC[:, :, :, :, 0] = np.diag(C0.flatten(1)).reshape(dim+1, dim+1,
dim+1, dim+1, order='F').copy()
zo13 = np.zeros((2, 1, 2, 1))
zo13[:,0,:,0] = np.mat([(1, 0), (0, 1)])
zo24 = np.zeros((1, 2, 1, 2))
zo24[0,:,0,:] = zo13[:, 0, :, 0]
pn13 = np.zeros((2, 1, 2, 1))
pn13[:,0,:,0] = np.mat([(1, -1), (-1, 1)]).reshape(2, 1, 2, 1,
order='F').copy()
pn24 = np.zeros((1, 2, 1, 2))
pn24[0, :, 0, :] = pn13[:, 0, :, 0]
for i in range(dim+1):
# Source conductance.
# First driver in the middle of row 1.
GG[dim/2, i, dim/2, i, 0] = G0
for j in range(dim):
# Horizontal segments.
node = 1 + j + i * dim
CC[j:j+2, i, j:j+2, i, node] = beta * zo13[:, 0, :, 0]
GG[j:j+2, i, j:j+2, i, node] = alpha * pn13[:, 0, :, 0]
# Vertical segments.
node = node + dim * ( dim + 1 )
CC[i, j:j+2, i, j:j+2, node] = beta * zo24[0, :, 0, :]
GG[i, j:j+2, i, j:j+2, node] = alpha * pn24[0, :, 0, :]
# Reshape for ease of use.
CC = CC.reshape((n*n, m+1), order='F').copy()
GG = GG.reshape((n*n, m+1), order='F').copy()
#
# Compute points the tradeoff curve, and the three sample points.
#
npts = 50
delays = np.linspace(50, 150, npts)
xdelays = [50, 100]
xnpts = len(xdelays)
areas = np.zeros(npts)
xareas = dict()
# Iterate over all points, and revisit specific points
for i in range(npts + xnpts):
# First pass, only gather minimal data from all cases.
if i < npts:
delay = delays[i]
print( ('Point {} of {} on the tradeoff curve ' \
+ '(Tmax = {})').format(i+1, npts, delay))
# Second pass, gather more data for specific cases,
# and make plots (later).
else:
xi = i - npts
delay = xdelays[xi]
print( ('Particular solution {} of {} ' \
+ '(Tmax = {})').format(xi+1, xnpts, delay))
#
# Construct and solve the convex model.
#
# Variables.
xt = cvx.Variable(m+1) # Element 1 of xt == 1 below.
G = cvx.Variable(n,n) # Symmetric constraint below.
C = cvx.Variable(n,n) # Symmetric constraint below.
# Objective.
obj = cvx.Minimize(cvx.sum_entries(C))
# Constraints.
constraints = [ xt[0] == 1,
G == G.T,
C == C.T,
G == cvx.reshape(GG*xt,n,n),
C == cvx.reshape(CC*xt,n,n),
delay * G - C == cvx.semidefinite(n),
0 <= xt[1:],
xt[1:] <= wmax,
]
#Solve problem
prob = cvx.Problem(obj, constraints)
prob.solve()
if prob.status != cvx.OPTIMAL:
raise Exception('CVXPY Error')
# Chop off the first element of x, which is
# constrainted to be 1
x = xt.value[1:]
# First pass, only gather minimal data from all cases.
if i < npts:
areas[i] = sum(x)
# Second pass, gather more data for specific cases,
# and make plots.
else:
xareas[xi] = sum(x)
#
# Print display sizes.
#
print 'Solution {}:'.format(xi+1)
print 'Vertical segments:'
print x[0:dim*(dim+1)].reshape(dim, dim+1, order='F').copy()
print 'Horizontal segments:'
print x[dim*(dim+1):].reshape(dim, dim+1, order='F').copy()
#
# Determine and plot the step responses.
#
A = -np.linalg.inv(C.value)*G.value
B = -A*np.ones((n, 1))
T = np.linspace(0, 500, 2000)
Y = simple_step(A, B, T[1], len(T))
indmax = -1
indmin = np.inf
for j in range(Y.shape[0]):
inds = np.amin(np.nonzero(Y[j, :] >= 0.5)[1])
if ( inds > indmax ):
indmax = inds
jmax = j
if ( inds < indmin ):
indmin = inds
jmin = j
tthres = T[indmax]
GinvC = np.linalg.solve(G.value, C.value)
tdom = max(np.linalg.eig(GinvC)[0])
elmore = np.amax(np.sum(GinvC.T, 0))
plt.figure(figsize=(8, 8))
plt.plot( T, np.asarray(Y[jmax,:]).flatten(), '-',
T, np.asarray(Y[jmin,:]).flatten() )
plt.plot( tdom * np.array([1, 1]), [0, 1], '--',
elmore * np.array([1, 1]), [0, 1], '--',
tthres * np.array([1, 1]), [0, 1], '--' )
plt.xlim([0, 500])
plt.ylim([0, 1])
plt.text(tdom, 0.92, 'd')
plt.text(elmore, 0.88, 'e')
plt.text(tthres, 0.96, 't')
plt.text( T[600], Y[jmax, 600], 'v{}'.format(jmax+1))
plt.text( T[600], Y[jmin, 600], 'v{}'.format(jmin+1))
plt.title(('Solution {} (Tmax={}), fastest ' \
+ 'and slowest step responses').format(xi+1, delay), fontsize=16)
plt.show()
#
# Plot the tradeoff curve.
#
plt.figure(figsize=(8, 8))
ind = np.isfinite(areas)
plt.plot(areas[ind], delays[ind])
plt.xlabel('Area')
plt.ylabel('Tdom')
plt.title('Area-delay tradeoff curve')
# Label the specific cases.
for k in range(xnpts):
plt.text(xareas[k][0, 0], xdelays[k], '({})'.format(k+1))
plt.show()