% L_cost_min.m    Find clc system with minimimum cost           20-8-15

% L_July3_Data1..6 files are merged into file L_July3 
load('D:\Power_Elect\Results\L\Opt_Data\L_Aug1');       % => l_data

% l_data(18,4,9,4,3,2,5,15)=(xd,km,tr,f,v12,fwd,pdes, elements).   
% (15) = xd,km,tr,f,v1,v2,phi,m1,m2,pin, pout,i1,ib2,vsm1,vsm2

char = 1e9*ones(18,4,9,4,4,18);     % h/w charact. xd, km, tr, f, kr(factor),                    
% elements (18) = xd, km, tr, f, kr, magmax(4), phs, eout, ein, ...
%               ecost, trcost, hscost, hwcost, totcost, phi_max.

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

hrs = 2000;
% Enter power useage profile:   runt(Vdc1/2(1:3), Pout(1:5), Fwd/Rev(1:2))
% CAUTION: AVOID P > 0.6 FOR 280/400 OR 400/280
%   Vdc1/2 (1:3) = 280/400, 400/400, 400/280;
%   Pout (1:5) = 0.2, 0.4, ... 1 * Pmax
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);    % energy out, kWhrs
        end
    end
end

% There are 5 H/W characteristic parameters: xd,km,tr,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 tr    
            for d = 1:4                 % range of freq
                exit = 0;               % set for pdes fail, => next hardware
% Reset the system maximum magnetics Irms and V.s : 
                magmax = zeros(1,4);                    % i1m,ib2m,vsm1m,vsm2m
                phi_max = 0;                            % to find max value
                for f = 1:3                             % v1v2 combinations
                    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:15) = l_data(a,b,c,d,f,h,g,:);	
                            if data(11) > 1           	% pdes achieved
% Record to largest value of phi. For each h/w set it should be > 0.36                                
                                if data(7) > phi_max
                                    phi_max = data(7);
                                end
% Find the largest (worst case) magnetics value for each system charact:
                                mag = data(12:15);     	% 2 Irms, 2 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) || (phi_max < 0.36)
                  	continue;
                end
% Trfmr Ivsp product = vsm1*i1rms + vsm3*ib2rms:
                ivsptr = magmax(3)*magmax(1) + magmax(4)*magmax(2); 
% Trfmr power loss is based on (IVsp)^1.17*kr*f^0.61
                ptr = 42e-9*ivsptr^1.17*data(4)^0.61;                    
% Power loss factor is the 5th hardware parameter:
             	for e = 1:4      	% range of 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:15) = l_data(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. Split Trfmr Pr between primary and secondary.
prtr = 0.58*ptr*((data(12)/magmax(1))^2 + (data(13)/magmax(2))^2)/2;  
% At optimal setting Pmag = 42% of mag component loss, is proportional to
%   Vs^2.8. Split Trfmr Pmag between primary and secondary.
pmagtr = 0.42*ptr*((data(14)/magmax(3))^2.8 + (data(15)/magmax(4))^2.8)/2;
% Calculate extra power loss (pex) from Trfmr.
pex = (prtr + pmagtr)*kr;                               % weighted extra p
ein = ein + (data(10) + pex)*runt(f,g,h);               % accumulate 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
              	trcost = 98e-16*ivsptr^2.1*data(4)^1.1/kr^1.8*kmc(b);
                hscost = varcost1*phs;
               	ecost = unitcost*esup;                	% energy cost
               	hwcost = fxdcost + trcost + hscost;   	% hardware cost
               	totcost = ecost + hwcost;
              	char(a,b,c,d,e,:) = [data(1:3) data(4)/1e3 kr magmax(1:2)...
                    magmax(3:4) phs eout/ein esup ...
                    ecost trcost hscost hwcost totcost phi_max];    % 18
                end                                  	% e, kr
            end                                         % d, freq
        end                                             % c, tr
    end                                                 % b, km
end                                                     % a, xd

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

for n = 1:10
    hwtc = char2(:,:,:,:,:,17);            	% 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, 18 elements
    result(:,n) = hw(:);
    char2(a,b,c,d,e,17) = 1e9;              % rule out entry in nxt srch
end





