% Digital circuit sizing for an inverter chain.
% (a figure is generated)
%
% This is an example taken directly from the paper:
%
%   Digital circuit optimization via geometrical programming
%   by Boyd, Kim, Patil, and Horowitz
%   Operations Research 53(6): 899-932, 2005.
%
% We consider a chain of N inverters driving a load capacitance CL.
% The problem is to find optimal scale factors for the inverter
% that minimize the sum of them (area), while obeying constraints
% on the maximum delay through the circuit, and minimum and maximum
% limits on scale factors. There are no limits on the total power.
% (For more details about the inverter chain see sec. 2.1.11 in the paper.)
%
%   minimize   sum(x)
%       s.t.   T_j <= Dmax          for j an output gate
%              T_j + d_i <= T_i     for j in FI(i)
%              x_min <= x <= x_max
%
% where variables are x and T.
% Here we use data structures and digital circuit models from the
% referenced paper.
%
% Almir Mutapcic 01/26/06
clear all; close all;

%********************************************************************
% problem data
%********************************************************************
N  = 8;      % number of inverters
CL = 20;     % capacitance load
Dmax = 20;   % maximum delay through the circuit
x_min = 1;   % minimum scale factor
x_max = 20;  % maximum scale factor

% primary inputs and primary outputs labels (start with N+1)
primary_inputs  = [N+1];
primary_outputs = [N+2];
M = N + length( primary_inputs ) + length( primary_outputs );

% fan-in cell array
FI{1} = [N+1];
for k = 2:N
  FI{k} = [k-1];
end
FI{N+2} = [N];

% fan-out cell array (compute it from the fan-in cell array)
FO = cell(M,1);
for gate = [1:N primary_outputs]
  preds = FI{gate};
  for k = 1:length(preds)
    FO{preds(k)}(end+1) = gate;
  end
end

% input and internal capacitance of gates and the driving resistance
Cin_norm  = ones(N,1);
Cint_norm = ones(N,1);
Rdrv_norm = ones(N,1);

% place extra capacitance before the input of the 5th inverter
Cin_norm(5) = 80;

% primary output has Cin capacitance (but has no Cload)
Cin_po = sparse(M,1);
Cin_po(primary_outputs) = CL;

% primary input has Cload capacitance (but has no Cin)
Cload_pi = sparse(M,1);
Cload_pi(primary_inputs) = 1;

%********************************************************************
% optimization
%********************************************************************
% optimization variables
gpvar x(N)                 % sizes
gpvar T(N)                 % arrival times

% input capacitance is an affine function of sizes
Cin  = Cin_norm.*x;
Cint = Cint_norm.*x;

% driving resistance is inversily proportional to sizes
R = Rdrv_norm./x;

% gate delay is the product of its driving resistance and load cap.
Cload = posynomial;
for gate = 1:N
  if ~ismember( FO{gate}, primary_outputs )
    Cload(gate) = sum( Cin(FO{gate}) );
  else
    Cload(gate) = Cin_po( FO{gate} );
  end
end
Cload = Cload';

% delay
D = 0.69*ones(N,1).*R.*( Cint + Cload );

% create timing constraints
timing_constr = [];
for gate = 1:N
  if ~ismember( FI{gate}, primary_inputs )
    for j = FI{gate}
      % enforce T_j + D_j <= T_i over all gates j that drive i
      timing_constr =  [timing_constr; D(gate) + T(j) <= T(gate)];
    end
  else
    % enforce D_i <= T_i for gates i connected to primary inputs
    timing_constr = [timing_constr; D(gate) <= T(gate)];
  end
end

% circuit delay is the max of arrival times for output gates
output_gates = [FI{primary_outputs}];
circuit_delay = max( T(output_gates) );

% collect all the constraints
constr = [timing_constr; circuit_delay <= Dmax;
          x_min*ones(N,1) <= x; x <= x_max*ones(N,1)];

% minimize the sum of scale factors subject to above constraints
[obj_value, solution, status] = gpsolve(sum(x), constr);
assign(solution);

% message about extra capacitance and result display
disp(' ')
disp(['Note: there is an extra capacitance between the 4th and 5th inverter'...
     ' in the chain.'])
fprintf(1,'\nOptimal scale factors are: \n'), x

% plot scale factors and maximum delay for inverter i
close all;
subplot(2,1,1); plot([1:N],T,'g--',[1:N],T,'bo');
ylabel('maximum delay T')
subplot(2,1,2); stem([1:N],x);
ylabel('scale factor x')
xlabel('inverter stage')
 
Iteration     primal obj.         gap        dual residual     previous step.
 
   1         5.04112e+00       5.70000e+01       2.49e+01               Inf
   2         6.26479e+00       2.34391e+01       8.19e-03       9.84878e-01
   3         5.65757e+00       2.04433e+01       4.32e-03       2.50000e-01
   4         4.94291e+00       1.60594e+01       1.27e-03       5.00000e-01
   5         4.35304e+00       1.17226e+01       6.21e-05       1.00000e+00
   6         1.71517e+00       5.85718e+00       1.77e-05       1.00000e+00
   7         6.20797e-01       2.93202e+00       1.24e-05       1.00000e+00
   8         2.47617e-01       1.47119e+00       7.91e-06       1.00000e+00
   9         1.17274e-01       7.39247e-01       4.29e-06       1.00000e+00
  10         4.42752e-02       3.71424e-01       8.69e-07       1.00000e+00
  11         4.54642e-03       1.86580e-01       2.45e-07       1.00000e+00
  12        -1.66740e-02       9.36691e-02       3.98e-08       1.00000e+00
 
Iteration     primal obj.         gap        dual residual     previous step.
 
   1         4.27906e+00       5.70000e+01       9.13e+01               Inf
   2         4.22971e+00       2.78393e+01       4.00e-02       1.00000e+00
   3         4.17149e+00       1.39168e+01       6.21e-03       1.00000e+00
   4         4.09986e+00       6.96255e+00       4.07e-03       1.00000e+00
   5         3.97608e+00       3.49228e+00       2.39e-03       1.00000e+00
   6         3.79703e+00       1.76392e+00       7.27e-04       1.00000e+00
   7         3.63606e+00       8.90384e-01       1.06e-04       1.00000e+00
   8         3.55343e+00       4.48155e-01       8.05e-06       1.00000e+00
   9         3.51320e+00       2.24976e-01       5.64e-07       1.00000e+00
  10         3.49332e+00       1.12746e-01       4.26e-08       1.00000e+00
  11         3.48341e+00       5.64429e-02       3.02e-09       1.00000e+00
  12         3.47846e+00       2.82395e-02       1.98e-10       1.00000e+00
  13         3.47598e+00       1.41243e-02       1.25e-11       1.00000e+00
  14         3.47474e+00       7.06333e-03       7.80e-13       1.00000e+00
  15         3.47412e+00       3.53195e-03       4.87e-14       1.00000e+00
  16         3.47381e+00       1.76605e-03       3.04e-15       1.00000e+00
  17         3.47366e+00       8.83042e-04       1.90e-16       1.00000e+00
  18         3.47358e+00       4.41526e-04       1.19e-17       1.00000e+00
  19         3.47354e+00       2.20764e-04       7.43e-19       1.00000e+00
  20         3.47352e+00       1.10382e-04       4.64e-20       1.00000e+00
  21         3.47351e+00       5.51912e-05       2.90e-21       1.00000e+00
  22         3.47351e+00       2.75956e-05       1.81e-22       1.00000e+00
  23         3.47351e+00       1.37978e-05       1.13e-23       1.00000e+00
  24         3.47351e+00       6.89891e-06       7.08e-25       1.00000e+00
  25         3.47351e+00       3.44945e-06       4.44e-26       1.00000e+00
  26         3.47351e+00       1.72473e-06       2.74e-27       1.00000e+00
  27         3.47351e+00       8.62363e-07       1.94e-28       1.00000e+00
  28         3.47351e+00       4.31182e-07       1.22e-29       1.00000e+00
  29         3.47351e+00       2.15591e-07       6.34e-30       1.00000e+00
  30         3.47351e+00       1.07795e-07       6.24e-30       1.00000e+00
  31         3.47351e+00       5.38977e-08       4.47e-30       1.00000e+00
  32         3.47351e+00       2.69489e-08       7.64e-30       1.00000e+00
  33         3.47351e+00       1.34744e-08       4.56e-30       1.00000e+00
  34         3.47351e+00       6.73721e-09       7.28e-30       1.00000e+00
Solved
Problem succesfully solved.
 
Note: there is an extra capacitance between the 4th and 5th inverter in the chain.

Optimal scale factors are: 

x =

    3.0461
    2.6553
    4.3323
   12.4395
    1.0000
    1.3185
    2.2358
    5.2220