% Digital circuit sizing example.
% (a figure is generated if the tradeoff flag is turned on)
%
% This is an example taken directly from the paper:
%
%   Digital Circuit Optimization via Geometric Programming
%   by Boyd, Kim, Patil, and Horowitz.
%   Operations Research 53(6): 899-932, 2005.
%
% Solves the problem of choosing gate scale factors x_i to give
% minimum ckt delay, subject to limits on the total area and power.
% Uses max gate arrival time T formulation that avoids evaluation
% of the delay over all paths in the circuit.
%
%   minimize   T_bar
%       s.t.   T_j <= T_bar      for j an output gate
%              T_j + d_i <= T_i  for j in FI(i)
%              P <= Pmax, A <= Amax
%              x >= 1
%
% where variables are x and T.
%
% We use the circuit topology presented in figure 1 (page 902),
% where we take gates 1, 3 and 6 to be inverters (INV),
% gates 2 and 7 to be three input NANDs (NAND3),
% and gates 4 and 5 to be two input NORs (NOR2).
%
% Almir Mutapcic 01/26/06
clear all; close all;
PLOT_TRADEOFF = 1; % to disable plotting set PLOT_TRADEOFF = 0;

%********************************************************************
% user specified data (specify problem constant and ckt topology)
%********************************************************************
m = 7;        % number of gates
Vdd = 5;      % supply voltage
Amax = 250;   % maximum area spec
Pmax = 15;    % maximum power spec

% gate specs
INV   = struct('Cin',3, 'Cint',3, 'Rdrv',0.48, 'A',3,  'Ileak',0.006);
NAND3 = struct('Cin',4, 'Cint',6, 'Rdrv',0.48, 'A',8,  'Ileak',0.007);
NOR2  = struct('Cin',5, 'Cint',6, 'Rdrv',0.48, 'A',10, 'Ileak',0.009);

gates([1 3 6]) = INV;
gates([2 7])   = NAND3;
gates([4 5])   = NOR2;

% primary inputs and primary outputs labels (start with m+1)
primary_inputs = [8 9 10];
primary_outputs = [11 12];
M = m + length( primary_inputs ) + length( primary_outputs );

% fan-in cell array
FI = cell(M,1);
FI{1} = [8];
FI{2} = [8 9 10];
FI{3} = [10];
FI{4} = [1 2];
FI{5} = [2 3];
FI{6} = [4];
FI{7} = [3 4 5];
FI{8} = [];
FI{9} = [];
FI{10} = [];
FI{11} = [6];
FI{12} = [7];

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

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

% activity frequency of gates and primary inputs
f_gates = 0.001*ones(m,1);
f_pi = sparse(M,1);
f_pi(primary_inputs) = 0.001*[10 10 10];

%********************************************************************
% derived problem data (computed from user inputs)
%********************************************************************
% fan-out cell array (compute it from the fan-in cell array)
FO = cell(M,1);
for gate = [1:m 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 driving resistance
Cin_norm  = [gates.Cin]';
Cint_norm = [gates.Cint]';
Rdrv_norm = [gates.Rdrv]';

% area specification for each gate with unit scaling
A_norm = [gates.A]';

% leakage current of gate i with unit scaling
Ileak_norm = [gates.Ileak]';

%********************************************************************
% optimization
%********************************************************************
% optimization variables
gpvar x(m)                 % scale factor
gpvar T(m)                 % 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:m
  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(m,1).*R.*( Cint + Cload );

% total area
area = A_norm'*x;

% total power calculation
Pdyn = Vdd^2*sum( f_pi(primary_inputs).*Cload_pi(primary_inputs) ) + ...
       Vdd^2*(f_gates'*(Cint + Cload));
Pstat = Vdd*Ileak_norm'*x;
power = Pdyn + Pstat;

% constraints
x_constr = ones(m,1) <= x; % all sizes greater than 1 (normalized)

% create timing constraints
timing_constr = [];
for gate = 1:m
  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

% objective is the upper bound on the overall delay
% and that is the max of arrival times for output gates
output_gates = [FI{primary_outputs}];
D = max( T(output_gates) );

% collect all the constraints
constr = [area <= Amax; power <= Pmax; timing_constr; x_constr];

% formulate the GP problem and solve it
[obj_value, solution, status] = gpsolve(D, constr);
assign(solution);

ckt_delay_plot = obj_value;
Pmax_plot = Pmax;

fprintf(1,'\nOptimal circuit delay for Pmax = %d and Amax = %d is %3.4f.\n', ...
        Pmax, Amax, obj_value)
disp('Optimal scale factors are: ')
x

%********************************************************************
% tradeoff curve code
%********************************************************************
if( PLOT_TRADEOFF )

% varying parameters for an optimal trade-off curve
N = 25;
Pmax = linspace(10,20,N);
Amax = [250];
min_delay = zeros(N,1);

% set the quiet flag (no solver reporting)
global QUIET; QUIET = 1;

for n = 1:N
  % add constraints that have varying parameters
  constr(1) = area  <= Amax(k);
  constr(2) = power <= Pmax(n);

  % solve the GP problem and compute the optimal volume
  [obj_value, solution, status] = gpsolve(D, constr);
  if strcmp(status,'Infeasible')
    error('Problem is infeasible.')
  end
  min_delay(n) = obj_value;
end

% enable solver reporting again
global QUIET; QUIET = 0;

% plot the tradeoff curve

plot(Pmax,min_delay); hold on
plot(Pmax_plot, ckt_delay_plot,'o'); hold off;
xlabel('Pmax'); ylabel('Dmin');
title(['Tradeoff curve for Amax = ' num2str(Amax)])

end
% end tradeoff curve code
 
Iteration     primal obj.         gap        dual residual     previous step.
 
   1         2.84042e+00       5.20000e+01       1.57e+02               Inf
   2         3.21190e+00       1.91274e+01       1.44e-01       9.74409e-01
   3         2.27880e+00       1.54051e+01       6.24e-02       3.59104e-01
   4         1.83068e+00       1.31067e+01       2.78e-02       3.50245e-01
   5         1.29195e+00       1.06673e+01       7.14e-03       5.24329e-01
   6         6.37228e-01       7.12541e+00       1.50e-03       5.89986e-01
   7         1.23548e-01       3.88893e+00       1.35e-04       7.90583e-01
   8        -4.09412e-02       2.24768e+00       5.14e-05       7.75006e-01
 
Iteration     primal obj.         gap        dual residual     previous step.
 
   1         1.79400e+02       5.20000e+01       1.09e+00               Inf
   2         1.29843e+02       4.96696e+01       1.03e+00       3.11976e-02
   3         9.38594e+01       4.85624e+01       9.87e-01       1.92534e-02
   4         6.93603e+01       4.73003e+01       9.32e-01       2.81792e-02
   5         4.87168e+01       4.54961e+01       8.48e-01       4.60102e-02
   6         3.10620e+01       4.27335e+01       7.15e-01       8.16028e-02
   7         9.33535e+00       3.45646e+01       4.04e-01       2.50000e-01
   8         8.24512e+00       3.00784e+01       1.06e-01       5.00000e-01
   9         8.01827e+00       2.59395e+01       2.76e-04       1.00000e+00
  10         5.07529e+00       1.30083e+01       1.25e-03       1.00000e+00
  11         3.70536e+00       6.56633e+00       5.49e-04       1.00000e+00
  12         2.98293e+00       3.31699e+00       2.30e-04       1.00000e+00
  13         2.61767e+00       1.67489e+00       7.76e-05       1.00000e+00
  14         2.43623e+00       8.44638e-01       2.84e-05       1.00000e+00
  15         2.34512e+00       4.26151e-01       1.46e-05       1.00000e+00
  16         2.29727e+00       2.15165e-01       5.99e-06       1.00000e+00
  17         2.27196e+00       1.08575e-01       1.56e-06       1.00000e+00
  18         2.25875e+00       5.46688e-02       2.49e-07       1.00000e+00
  19         2.25203e+00       2.74571e-02       2.51e-08       1.00000e+00
  20         2.24864e+00       1.37643e-02       1.82e-09       1.00000e+00
  21         2.24693e+00       6.89201e-03       1.18e-10       1.00000e+00
  22         2.24607e+00       3.44859e-03       7.44e-12       1.00000e+00
  23         2.24564e+00       1.72496e-03       4.67e-13       1.00000e+00
  24         2.24542e+00       8.62646e-04       2.92e-14       1.00000e+00
  25         2.24532e+00       4.31365e-04       1.82e-15       1.00000e+00
  26         2.24526e+00       2.15693e-04       1.14e-16       1.00000e+00
  27         2.24524e+00       1.07849e-04       7.12e-18       1.00000e+00
  28         2.24522e+00       5.39252e-05       4.45e-19       1.00000e+00
  29         2.24522e+00       2.69628e-05       2.78e-20       1.00000e+00
  30         2.24521e+00       1.34814e-05       1.74e-21       1.00000e+00
  31         2.24521e+00       6.74072e-06       1.09e-22       1.00000e+00
  32         2.24521e+00       3.37036e-06       6.79e-24       1.00000e+00
  33         2.24521e+00       1.68518e-06       4.25e-25       1.00000e+00
  34         2.24521e+00       8.42592e-07       2.65e-26       1.00000e+00
  35         2.24521e+00       4.21296e-07       1.66e-27       1.00000e+00
  36         2.24521e+00       2.10648e-07       1.10e-28       1.00000e+00
  37         2.24521e+00       1.05324e-07       1.06e-29       1.00000e+00
  38         2.24521e+00       5.26620e-08       2.07e-30       1.00000e+00
  39         2.24521e+00       2.63310e-08       2.94e-30       1.00000e+00
  40         2.24521e+00       1.31655e-08       7.37e-30       1.00000e+00
  41         2.24521e+00       6.58275e-09       1.33e-30       1.00000e+00
Solved
Problem succesfully solved.

Optimal circuit delay for Pmax = 15 and Amax = 250 is 9.4424.
Optimal scale factors are: 

x =

    3.2339
   13.6461
    3.9632
    3.2337
    2.2241
    1.3157
    2.1737

Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.
Problem succesfully solved.