% lcl_cost_min.m    Find lcl system with minimimum cost         29-9-15

% LCL_Jun30_Data1..6 plus 200 kHz files are merged into file LCL_Sept24 
load('D:\Power_Elect\Results\LCL\Opt_Data\LCL_Sept24');      % => lcldata

% lcldata(18,4,9,5,3,2,5,18)=(xd,km,trr,f,v12,fwd,pdes, 18 elements).   
% (18) = xd,km,trr,f,v1,v2,phi,m1,m2,pin, pout,i1,i3,ib2,ic1,vsl1,vsm1,vsm2

char = 1e9*ones(18,4,9,5,4,22);  % h/w charact. xd, km, trr, f, kr(factor)
% elements (22) = xd, km, trr, freq, kr, magmax(7), phs, eout, ein, ...
%             	ecost, lcost, trcost, hscost, capcost, hwcost, totcost.

unitcost = 0.29;             % cost of 1 kWHr of electricity
% Specify converter costs, assuming approximately 1000 units made. 
% Fixed (run power independent) and variable (varies with run pwr) costs: 
fxdcost = 150;              % fpga,driver,p/s,adc,pcb,transis,misc,dsgn
                            %   5,10,10,5,20,30,20,20 + 30 labour 
varcost1 = 0.1;           	% variable cost coeff 1: HS, $/phs (W)
% Assumes heatsink cost is prop. to its maximum power dissipation phs.
% For magnetic components, L1 and trfmr, assume:  
%   Cost is prop. to (IVsp)^2.1*f^1.1*kr^-1.8*[kmc]  
%   Power Loss is prop. to (IVsp)^1.17*f^0.61*kr.
%  	These are based on trfmr having 96k V.us.Arms @ 50 kHz, 21 W loss, 
%   	cost $41.4. Inductor having 40k V.us.Arms 11W loss, cost $21.2.
kmc = [0.85 0.9 0.95 1];   	% volume/cost factor for km
kmc = [1 1 1 1];   	% volume/cost factor for km

hrs = 2000;
% Power useage profile:   runt(Vdc1/2(1:3), Pout(1:5), Fwd/Rev(1:2))
%   Vdc1/2 (1:3) = 280/400, 400/400, 400/280; Pout (1:5) = 0.2 .. 1*Pmax
%   Enter relative hours of duty
runt = zeros(3,5,2);                 	% rel hrs = f(volts,power,fwd/rev)
runt(2,5,1) = 85;                       % 400/400, P, fwd
runt(2,5,2) = 10;                   	% 400/400, P, rev
% runt(1,1,1) = 2;                        % 280/400, 0.2P, fwd
runt(3,1,1) = 5;                        % 400/280, 0.2P, fwd
% runt(3,1,2) = 5;                        % 400/280, 0.2P, rev
% runt(1,4,2) = 6;                        % 280/400, 0.8P, rev
runt = hrs*runt/sum(sum(sum(runt,3)));       	% distribute hrs

eout = 0;
for a = 1:3
    for b = 1:5
        for c = 1:2
            eout = eout + b*0.8*runt(a,b,c);    % total energy out of
        end                                     %   converter, kWhrs
    end
end

% There are 5 H/W characteristic parameters: xd,km,trr,f. (The 5th parameter,  
%   the kr factor, will be used seperately)
for a = 1:18                            % range of xd
    for b = 1:4                         % range of km                           
        for c = 1:9                     % range of trr    
            for d = 1:4                 % range of freq
% Find worst case reactive (magnetic + capacitive) component required ratings.
                exit = 0;               % set for pdes fail, => next hardware
% Reset the system maximum magnetics Irms and V.s : 
                magmax = zeros(1,7);    % i1m,i3m,ib2m,ic1m,vsl1m,vsm1m,vsm2m
                for f = 1:3                             % v2 = 280,340,400
                    if exit > 0;
                        break;
                    end
                    for g = 1:5                     	% pdes 20% steps
                        if exit > 0;
                            break;
                        end
                        for h = 1:2                  	% fwd, rev
                            if runt(f,g,h) == 0     	% no pwr reqd, do nxt 
                                continue;             	%   v,p,dir setting 
                            end                                     
                            data(1:18) = lcldata(a,b,c,d,f,h,g,:);	
                            if data(11) > 1           	% pdes achieved
% Find the largest (worst case) magnetics value for each system charact:
                                mag = data(12:18);     	% 4 Irms, 3 V.us
                                magmax(mag > magmax) = mag(mag > magmax);
                            else                      	% pdes not met, go 
                                exit = 1; break;      	% 	to next h/w
                            end
                        end
                    end
                end
              	if exit > 0;
                  	continue;
                end
% Calculate C reactive size :
                capv = magmax(4)^2*data(1)/data(4);     % I^2.Xd/f
                ivspl1 = magmax(5)*magmax(1);         	% L1 Ivsp = vsl1*i1rms:                                                             
% Trfmr Ivsp product = vsm1*i3rms + vsm2*ib2rms:
                ivsptr = magmax(6)*magmax(2) + magmax(7)*magmax(3); 
% Calculate L1 and Trfmr power loss for maximum load and reduction factor = 1
                ptr = 42e-9*ivsptr^1.17*data(4)^0.61;                    
                pl1 = 62e-9*ivspl1^1.17*data(4)^0.61;	
% Power loss factor is the 5th hardware parameter:
             	for e = 1:4      	% range for pwr loss factor, 
                    if exit > 0;
                        break;
                    end
                    kr = 2^(1 - e);	% pwr reduction factor 1/1,2,4,8   
% Add additional losses (from magnetic components) and accumulate energy.
                    exit = 0;     	% set for pdes fail, => next hardware
                    ein = 0;      	% Reset the total input energy 
                    esup = 0;     	% Reset the total energy supplied
                    phs = 0;       	% Reset the system maximum heatsink power 
                    for f = 1:3                       	% v2 = 280,340,400
                        if exit > 0;
                            break;
                        end
                        for g = 1:5                   	% pdes 20% steps
                            if exit > 0;
                                break;
                            end
                            for h = 1:2              	% fwd, rev
                                if runt(f,g,h) == 0   	% no pwr reqd, do nxt 
                                    continue;          	% 	v,p,dir setting
                                end       
                                data(1:18) = lcldata(a,b,c,d,f,h,g,:);	
                                if data(11) > 1      	% pdes achieved
%------------------------------------------------------------------------
% At optimal setting Pr = 58% of mag component loss. Pr is proportional to
%   I^2 at other settings. In trfmr split Pr between primary and secondary.
prl1 = 0.58*pl1*(data(12)/magmax(1))^2;                 % Pr prop. to I^2
prtr = 0.58*ptr*((data(13)/magmax(2))^2 + (data(14)/magmax(3))^2)/2;  
% At optimal setting Pmag = 42% of mag component loss, is proportional to
%   Vs^2.8. Split Trfmr Pmag between primary and secondary.
pmagl1 = 0.42*pl1*(data(16)/magmax(5))^2.8;             % Pmag prop to Vs^2.8
pmagtr = 0.42*ptr*((data(17)/magmax(6))^2.8 + (data(18)/magmax(7))^2.8)/2;
% Calculate extra power loss (pex) from mag components L1 and Trfmr.
pex = (prl1 + prtr + pmagl1 + pmagtr)*kr;               % weighted extra p
ein = ein + (data(10) + pex)*runt(f,g,h);             	% accum input energy
% Calculate the net energy input from the supply:
if (h == 2)                                          	% for reverse
    esup = esup - data(11)*runt(f,g,h);               	% pout goes back to mains
else
    esup = esup + (pex + data(10))*runt(f,g,h);        	% extra + pin taken from mains
end
loss = data(10) - data(11) + pex;                       % pin - pout + extra
if loss > phs                                           % find maximum dissip
  	phs = loss;                                         %   of heatsink
end
%------------------------------------------------------------------------                    
                                else                    % pdes not met, do next                                     
                                    exit = 1; break;   	%	hardware
                                end
                            end                         % h, direction
                        end                          	% g, pdes
                    end                               	% f, v12
              	if exit > 0;
                  	continue;
                end
                esup = esup/1e3;  ein = ein/1e3;       	% now kWHr
                lcost = 31e-15*ivspl1^2.1*data(4)^1.1/kr^1.8;
                trcost = 98e-16*ivsptr^2.1*data(4)^1.1/kr^1.8*kmc(b);
                hscost = varcost1*phs;
                capcost = 2.5 + 2*capv;
               	ecost = unitcost*esup;                  % energy cost
               	hwcost = fxdcost + lcost + trcost + hscost + capcost;  	% hardware cost
               	totcost = ecost + hwcost;
              	char(a,b,c,d,e,:) = [data(1:3) data(4)/1e3 kr  ...
                    magmax(1:4) magmax(5:7) phs eout/ein esup ...
                    ecost lcost trcost hscost capcost hwcost totcost];    % 22
                end                                  	% e, kr
            end                                         % d, freq
        end                                             % c, trr
    end                                                 % b, km
end                                                     % a, xd

char2 = char;                                       	% make working copy
result = zeros(22,10);

for n = 1:10
    hwtc = char2(:,:,:,:,:,22);            	% each hw total cost, dim = 5 
    [mincost,ind] = min(hwtc(:));          	% linear indice
    [a,b,c,d,e] = ind2sub(size(hwtc),ind);	% array addr of best h/w       
    hw = char2(a,b,c,d,e,:);              	% best hardware, 22 elements
    result(:,n) = hw(:);
    char2(a,b,c,d,e,22) = 1e9;              % rule out entry in nxt srch
end

