function litzopt()
% Returns transformer or inductor stranding designs of minimal cost and loss
%
% LITZOPT returns a table of AWG wire sizes and
% strand number and plots these designs relative
% to the optimal AWG 44 design. LITZOPT is for
% use only with piecewise-linear current waveforms.
% The user runs LITZOPT to create a data file without
% transformer chracteristic data. The user inputs
% the transformer data % into that file includingr
% the output (magnetic field value) of a finite
% element analysis simulation (e.g. Maxwell). The
% user then reruns LITZOPT to optimize
% the selection of diameter and number of litz
% strands considering cost and loss.
% Version 1.7 updated by Jennifer Pollock
% Updated functions testnan
% Updated buildbletable to handle 10 windings
% Version 1.6.5 updated by Jennifer Pollock
% Removed bugs so litzopt runs on matlab version 7.2
% Improved commenting to separate functions and explain some
% functionality in foreachecitation
% Version 1.6.4
% Jennifer Pollock
% Fixed the "max" index bugs that result when the code is run
% in Matlab Version v.7
%Version 1.6.3
%Jennifer Pollock
% 7/8/04 fixed testnan, testnandia, and testnanloss functions to scan the entire matrix for nan
%Version 1.6.2
% 9/3/03 fixed bugs in data entry no legs gapped option.
% 9/3/03 organized data entry in file created to 2d option
%Version 1.6.1
% 8/11/03 added windinght1d function to reduce 1d data entry
% 8/11/03 fixed fea function to generate a new data file - works now
%Version 1.6
% 8/2/03 added the ability plot an orignal design on the optimal design frontier
% and print cost and loss for original design to 1d, 2d, fea analysis
% 8/2/03 changed 1d data file to take in the bobbin area, 'areabb', set breadth of bobbin
% window equal to breadth of core window, set bobbin height, h = areabb/bw and set 'a' equal to
% the bobbin height divided by the number of windings - a = h/numwind
% 7/30/03 added function roundtoAWGrange that takes any AWG vector
% and rounds the dia matrix to values in AWG vector
% 7/17/03 jdp switched gap = b or c in B2avgcalc1d
%
% Version 1.5
% Jennifer Pollock
% updated the local version with slope matching optimization 6/22/03
%
% Version 1.4.1 5/1/03
% Charlie Sullivan
% Corrected comments on 1.4 and added semicolons where missing
% Corrected gap interpretation to do all legs gapped for 'A...' Removed 2D outer legs only, since that isn't supported for 2D.
%
% Version 1.4 4/24/03
% Charlie Sullivan
% Note: gaplen is length of each gap, not total.
% Updates: looked for best of 1.3, 1.29. Seems to be just cosmetic changes.
% Added space between matrix elements to avoid annoying error messages
% line 1396 in B2avgcalc2d, made numwind plus gap = numwind +2 where it was 1, and other changes in ungap
% Made centergap at left of window.
% removed length(yzero...) and replaced with correct dimension, named n1yzp. Line 1540 and the half dozen after that.
% fixed gap excitation magnitude in foreachexcitation
% note: does rectfldcalc for unexcited windings, with zero as amplitude--waste of time for larger number of windings--fixed--see below
% fixed B^2 = Bx^2 + by^2 not (bx +by)^2
% fixed Hx, Hy calc to accumulate over different windings, not throw out old data.
% Fixed imageclac to trim zeros out of list "excwind" before sending to Rectfldcalc
% Fixed rectfld calc to use an modified excwind that is the same size as the aexc, bexc, etc. variables--one element for each excited winding
% or image of an excited winding, instead of one element per orginal winding. This corrects the calc and eliminates an internal loop, thus
% speeding it up. Also and did some other little stuff to speed up, like taking 2*pi multiplication out of inner looop.
% Tested on 30:30 design from original cost/loss paper
% Note: imaginary answers typically result from diagonal elements too large from problem wiht field calc.
% 1.29 updates: at least some dialog box text improved
% 1.2 Updates; Simualtions now called internal instead of analytic,
% Typos for winex to windex in 1d code fixed
% Units specified in dialogue boxes
% All legs, center leg, and outer legs options used to describe
% gap positions in dialogue boxes
% Both outer legs case corrected (excitation vector had too many elements)
% By Dr. Charles R. Sullivan and Tarek Abdallah and Jennifer Pollock
% Dartmouth College
% www.thayer.darmouth.edu/inductor
% Please forward questions and comments to charles.sullivan@dartmouth.edu
% The winding loss estimation algorithm employed in
% this program is introduced in the paper "Winding
% Loss Calculation with Multiple Windings, Arbitrary
% Waveforms and Two-Dimensional Field Geometry"
% by Dr. Charles R. Sullivan from % IEEE Industry
% Applications Society Annual Meeting, Oct. 1999, Phoenix,
% pp. 2093-2099. (This paper is superceded by the 2001
% Transactions on Power Electronics article, but the
% original is still available if you really want
% it. )
%
%% The stranding cost estimation algorithm is from
% "Cost-Constrained Selection of Strand Wire and
% Number in a Litz-Wire Transformer Winding." by
% Dr. Charles % R. Sullivan IEEE Transactions
% on Power Electronics, 16 (2), March 2001,
% pp. 281-8. from "Proceedings of the 1998 IEEE
% Industry Applications Society
% % Annual Meeting, 1998, pp. 907-911.
% List Of Functions:
%
%
% % New: This file creates a new transformer data
% file with appropriate matrix and vector dimensions,
% but no data, to be edited by the user.
%
%
% Twodimsimulation: A function that is called when a two-dimensional field simulation is
% to be performed. The function prompts the user for data to be
% stored in a data file. The function executes a succession of commands
% to ask the user for various transformer data in a GUI format. The
% function concludes by closing the data file and saving it under the
% filename and directory of the user's choice. The data file is edited
% by the user and again called when a two-dimensional field simulation
% is to be performed to allow for subsequent strand optimization.
%
% Onedimsimulation: A function that is called when a one-dimensional field simulation is
% to be performed. The function prompts the user for data to be
% stored in a data file. The function executes a succession of
% commands to ask the user for various transformer data in a GUI
% format. The function concludes by closing the data file and
% saving it under the filename and directory of the user's choice.
% The data file is edited by the user and again called when a
% one-dimensional field simulation is to be performed to allow
% for subsequent strand optimization.
%
% FEA: A function that is called when the user requests field data to be externally
% computed by finite element analysis and inputted into a Litzopt data file.
% The function prompts the user for data to be stored in a data file. The function executes
% a succession of commands to ask the user for various transformer data in a GUI
% format. The function concludes by closing the data file and saving it under
% the filename and directory of the user's choice. The data
% file is edited by the user to include the results of a fea-based field
% simulation to allow for subsequent strand optimization.
%
%
% Existing: This file optimizes the stranding design
% given a file edited by the user
%
% averagebsquare: Generates average value of the squared field
% with respect to each of the windings due to the
% current in all the windings. One input is a three-dimensional
% array containing square of the sum of squared
% fields averaged over winding 'k' (for terms in page 'k')
% due to exited windings 'i' and 'j' (the output of a
% finite element analysis or analytic field approximation) on the
% main diagonal and upper tringle. This input has zero
% terms in lower traingle. Additional input is 'N' (turns number
% in each winding). The output is also a three-dimensional
% array of products Bi*Bj at each location (i,j) multiplied by
% the prodict of turns number Ni*Nj at location (i,j) on each page.
%
%
% %Costnloss: Returns vectors of costs, losses, and stranding
% for each of the calculated designs where each
% row designates a different design. Inputs are
% 'Irms' (root mean squared current), 'lt' (turn
% length), 'N' (number of turns), 'awg' (a vector
% of strand diameters for each of the optimized
% designs), 'rhoc' (the resistivity of copper at
% operating temperature), and 'kl' (a parameter
% defined by time-rate of change of currents and
% the modified dynamic loss matrix). The outputs
% are an array called 'cost' (Cost function of litz
% wire is assumed to be proportional to
% '(1+(k/(DC^4))+(kk))*n*lt*N' where 'k' is 1.1e-26,
% 'kk' is 2e-9, 'n' is the number of strands in the
% design in question 'lt' is the length of a single
% turn, and 'N' is the number of turns. where, dc
% is strand diameter, n is strand number, lt is turn
% length, and N is turn number), 'actual_loss'
% (the total winding power lost in each design and is
% equal to 'Irms^2*lt*N*rhoc*Fr/n/(pi*dc^2/4)' where
% all parameters are as defined above and rhoc is the
% resistivity of copper at operating temperature),
% and 'n' is defined to be the number of strands
% implemented in each design.
%
% %dCostm: Approximates the derivative of cost function given
% strand size and the same cost function is used
% (see Costnloss)
%
% Numstrands: Solves for the number of strands in the optimal
% design given strand size and ac to dc resistance
% ratio (Fr)
%
% Awgtomet: Converts from American Wire Gauge strand size measure
% to meters
%
% COSTm: Calculates cost basis proportional to the additional
% cost per unit mass of a given strand size. Cost is
% considered with the function given above (see Costnloss)
%
%
%
% B2avgcalc1d: Solves for magnitude of magnetic field in one spatial dimension.
% Inputs are breadth of winding window, gap locations,uo, and
% number of windings
%
%
%
%
%
% B2avgcalc2d: Generated two-dimensional field solution and output is area-integral
% of squared field, givens are transformer geometry such as, gap
% locations and size, winding locations and size, number of requested
% images and winding divisions. Calls function foreachexcitation for
% specific excitation
%
%
%
% foreachexcitation: changes winding excitations in sequence to obtain all field
% values described in the SFD method. inputs are winding positions
% and locations in window, gap length and number, breadth and height
% of the winding window, gap locations, winding divisions for field
% computation, vector describing excitation and number
% of images to be computed
%
% imagecalc: computes all data needed to use method of images about high perbeablilty
% core, the inputs are dimension and location of each of the windings in
% the window, breadth and height of the window, number of divisions for
% windings, excited winding vector, and number of images to be computed.
% The entire system of imaged excited windings are generated.:
%
% rectfldcalc: takes in the imaged system of excited windings and computes the
% area-average squared magnetic field over the cross-section of a
% single transformer winding. The inputs are the location and
% dimensions of the imaged system, the location and dimensions
% of the winding being averaged over, breadth and height of the
% winding window, and the number of divisions of the winding.
% The single area-avergae field value is the output.
% ungap: to place Ribbons with same length of gaps and magnitude
% equal and opposite to winding current; then ungap the core.
%
% PWLcurrentplot: takes in 'dt' (vector of time segment lengths for current waveform),
% 'I' (matrix of current current values), and 'numwind' to create current
% waveform plots
%
% Imscurrent: computes mean-square current Ims and root mean squared current Irms
% from I, dt, numtime, numwind
%
% find_k_ell: calculate di/dt for each winding and approximate loss by SFD method
% from numwind,numtime,I,Ims,dt,N,len,rhoc,B2avg and returns 'kl'
%
% find_optimal_designs: finds the designs for wire sized in vector awg, and also the
% set of full-bobbin designs (all the "2" ones)
% numwind, kl, Irms, len, N, rhoc are all pretty obvious.
% userfp is the maximum packing factor relative to perfect
% square packing of cylinders that the user specifies; used for finding
% full-bobbin designs. Takes in: awg,numwind,kl,Irms,len,N,rhoc,bb,a,insulbuild,userfp
% Returns:[totcost,totloss,lossW,ncorrect,fp,totcost2,totloss2,loss2,n2,totcost3,
% totloss3,loss3,n3, awground, n2round,totcost2round, totloss2round, dia
%
% find_optimal_designs_curvedata: Takes in awgfine,numwind,kl,Irms,len,N,rhoc,bb,a,insulbuild,userfp.
% Does the same thing as find_optimal_designs, but with awgfine to create
% data to plot. Only returns totcostfine,totlossfine,totcost2fine,totloss2fine
% where totcostfine and totlossfine are the hypothetical optimal design curve
% and totcost2fine and totloss2fine are the buikdabke design curve
%
% designlabels: creates designlabels for buildable design table
%
% plotdesignL: plots the optimal design frontier
%
% Xsectionplot: plots cross-section of winding window to confirm location of windings in the window
%
% findtable: This function is called by find_optimal_designs and find_optimal_design_curvedata.
% Finds the lookup table
%
% numderiv: This function calculates the numerical deriv by taking x,y or cost, loss and
% returning the derivative and it's location.
%
% fprime: This function looks up the derivtive in the look up table givena cost and returns the derivative
%
% fprimeinverse: This function looks up the normalized cost in the look up table given a derivative
%
% dlookup: This function looks up the diameter in the look up table given a cost
%
% dlookupinverse: %this function looks up the normalized cost cooresponding to a diameter(awg) in the look up table.
%
% ff: This the forward function that given a noralized cost, it looksup the normalized loss
%
% stranding: This function takes in Irms,lt,kl,N,rhoc,d2,lossW2 and returns 'n', the number of strands
% This function calls numstrand
%
% totaldesignsAWG: This function is called in find_optimal_designs. It takes in dia, awg, numwind
% and returns rounded designs 'diaround' matrix of round diameters.
% this function will determine all the rounded designs
% the criterion are for numbers with a remainer < 0.2, it is rounded down.
% If the remainder falls between 0.2 and 0.8, it is rounded both ways
% if the remainder is >0.8 the number id rounded up
%
% optpackfactor: This function will determine the packing factor of optimal the optimal
% design considering the insulation build. It returns the matirx 'fp' and
% takes in awg, N, n, bb, a, insulbuild
%
% buldabledesigns: This function will determine the a buildable design considering the insulation build
% It returns the matrix 'n2' of buildable design stranding
% It takes in awg, a, n, fp, N, userfp, bb, insulbuild
%
% totaldesigns: This function will determine all the rounded designs for rounding
% the stranding. The criterion are for numbers with a remainer < 0.2,
% it is rounded down. If the remainder falls between 0.2 and 0.8, it
% is rounded both ways. if the remainder is >0.8 the number is rounded up
% It returns awground, n2round and takes in n2, awg, numwind
%
% singlebuild: This function deterimes wire diameter with single build insulation
%
% heavybuild: This function determines wire diameter with heavy build insulation
%
% smartround: This function allows a number or vector to be rounded to a specified
% number of significant digits
%
% costnloss2: This is a modified version of Costnloss to determine the cost and
% losses for buildable designs the method remains much the same,
% accept the n2 matrix is passed in to the function and it returns
% just the cost2 and actual_loss2 vectors.
%
% fullbobbindesigns: this function will determine the full-bobbin designs using
% the insulation build
%
% optdesigntable: This function generates html table for optimal design given
% optdesignlabels, totcost, totloss, numwind, n, dia, fp
%
% buildabletable: This function generates Html table for buildable design
% uses the smartround function to round cost and loss columns
% It takes in: bbdesignlabels, totcost2round, totloss2round,
% numwind,n2round,awground
%
% rounddesigns: This function determines how to round designs and number of
% designs to create. It returns nrounded, awground and takes in awg,n2, numwind
%
% testnan: This function prunes matrix that contains NaN that result from
% trying to extrapolate from the lookup table/ It returns inputvector, kcostmatrix, klossmatrix
% and takes in inputvector, kcostmatrix, klossmatrix
%
% testnanloss: %This function prunes matrix that contains NaN that result
% from trying to extrapolate from the lookup table. It returns:
% inputvector, costW, kcostmatrix, klossmatrix and takes in:
% inputvector, costW, kcostmatrix, klossmatrix)
%
% testnandia: %This function prunes matrix that contains NaN that result from
% trying to extrapolate from the lookup table. It takes in input vector,
% costW, lossW, kcostmatrix, klossmatrix. It returns: inputvector, costW,
% lossW, kcostmatrix, klossmatrix
%
% roundtoAWGrange: this function rounded the matrix of wire diameters to the nearest
% value in awg vector passed in. It returns the matrix of rounded
% wire sizes.
%
% originaldesign: this function determines the cost and loss for the orignal design
% if data is entered into the fields 'orig_numstrands' and 'orig_wiresize'.
% It returns the total relative cost and loss for original design
%
% windinght1d: this function calculates h, ht of core window, hb, ht of bobbin window
% bb, breadth of bobbin window and a, ht of each winding using N-number of turns
% Irms - RMS currnet, bw - core breadth and areabb-area of bobbin window
%
%
%%%%%%%%%%%%%Start of Main Program%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
answer=questdlg('Select Project File:','Litzopt','New Project: Internal Field Simulation','New Project: External Field Simulation','Existing Project','New Project: Internal Field Simulation');
% Decide if Existing file or New,
if answer(1)=='E'
existing
elseif answer(14)=='I'
new
else
fea
end
%%%%%%%%%%%%%End of Main Program%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%start of new %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function new()
%this function setup a new data file for either a 1D or 2D solution
answer=questdlg('Geometry Considered:','Litzopt','1-D Field Geometery Considered','2-D Field Geometery Considered','2-D Field Geometery Considered');
if answer(1)=='1'
onedimsimulation
else
twodimsimulation
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end of new%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% twodsimulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function is run when a new data file is to be created and the
% user selects a 2 dimensional simulation. For the analytic sol'n, the input parameters no longer
% include winding volumes and simulated field values as they are calculated
% in another function. For the 2d case, extra user inputs ars gaplen, xzero, yzero (winding loc's),
% xdiv,ydiv (core divisions)
function twodimsimulation()
prompt={'Number of windings:','Number of time segments for piecewise-linear current waveforms:'};
answer=inputdlg(prompt);
vect=str2double(answer);
numtime=vect(2);
numwind=vect(1);
uo=4e-7*pi;
answer1=questdlg('Are there any gaps in the core? ','Choose','No legs gapped','Center leg gapped','All legs gapped','No legs gapped');
if answer1(1)=='C'
gaploc='Center';
elseif answer1(1) == 'A'
gaploc ='All';
else
gaploc ='No Gap';
gaplen = 0;
end
answer2=questdlg('What build is the wire insulation?','Choose','Single','Heavy','Single');
if answer2(1)=='S' %this is single build insulation
insulbuild='s';
else %this is heavy insulation
insulbuild = 'h';
end
% check to see if '.m' follows filename
[filename,pathname]=uiputfile('Untitled', 'Save new file as:');
both= and (filename(length(filename)-1)=='.',filename(length(filename))=='m');
% if '.m' does not follow filename, place it there
if both==0; % if false statements
filename=[filename '.m'];
end
% store original pathname and set path to new pathname
origdirect=cd;
cd(pathname);
% create new file in the desired path
[fid,errormsg]=fopen(filename,'w'); % %#ok<*MSNU>
% begin writing file
startstr='%%For an example and help, go to www.thayer.dartmouth.edu/inductor/litzopt, or refer to winding data for use with litzopt. \n \n';
startstr=[startstr '%%Please replace the question marks with the pertinent data and leave all else alone....\n \n'];
startstr=[startstr '%%There are ' num2str(numwind) ' windings in your transformer. \n \n'];
startstr=[startstr '%%Resistivity (1.77e-8 is for room temperature copper) \n '];
startstr=[startstr ' rhoc = 1.77e-8; %% [ohm-meters]\n \n'];
startstr=[startstr '%%Perform optimization on designs with wire sizes ranging from American Wire Gauge strand size (32 to 48 is the default) \n'];
startstr=[startstr ' awg = [32:2:48]; %% \n \n'];
startstr=[startstr '%%Number of images of the Winding Window Geometry to be computed in the horizontal direction (5 images is the default) \n'];
startstr=[startstr ' Ximage = 5; %% \n \n'];
startstr=[startstr '%%Number of images of the Winding Window Geometry to be computed in the vertical direction (5 images is the default) \n'];
startstr=[startstr ' Yimage = 5; %% \n \n'];
startstr=[startstr '%%Number of horizontal divisions of a winding (20 divisions is the default) \n'];
startstr=[startstr ' xdiv = 20; %% \n \n'];
startstr=[startstr '%%Number of vertical divisions of a winding (20 divisions is the default) \n'];
startstr=[startstr ' ydiv = 20; %% \n \n'];
startstr=[startstr '%%Breadth of the Core Window in meters \n'];
startstr=[startstr ' bw = ?? ; %% \n \n'];
startstr=[startstr '%%Height of the Core Window in meters \n'];
startstr=[startstr ' h = ?? ; %% \n \n'];
startstr=[startstr '%%Breadth of the Bobbin Window in meters \n'];
startstr=[startstr ' bb = ?? ; %% \n \n'];
startstr=[startstr '%%Height of the Bobbin Window in meters \n'];
startstr=[startstr ' hb = ?? ; %% \n \n'];
startstr=[startstr '%%Length of the gap in meters \n'];
startstr=[startstr ' gaplen = ?? ; %% \n \n'];
startstr=[startstr '%%Maximum Achieveable Packing Factor \n'];
startstr=[startstr ' userfp = ?? ; %% \n \n'];
fprintf(fid,startstr);
%%% Length, create vector of winding lengths
istr=['%%Average turn length for each winding... ' num2str(numwind) ' elements are required in this vector. \n'] ;
istr=[istr ' len =[ '];
%%%% loop through windings and designate vector element for each winding length
for i=1:numwind
istr=[istr ' ?? '];
end
% close vector and write it to file
istr=[istr ']; %%[meters] \n \n'];
fprintf(fid,istr);
%%% Time Duration %%%
% create vector of time durations for each segment
dtstr=['%%Duration of each time segment... ' num2str(numtime) ' elements are required in this vector. \n'];
dtstr=[dtstr ' dt =[ '];
% loop through number of time segments and designate a vector element for additional
% time segment duration
for j=1:numtime
dtstr=[dtstr ' ?? ']; %#ok<*AGROW>
end
% close vector and write it to file
dtstr=[dtstr ']; %%[seconds] \n \n'];
fprintf(fid,dtstr);
%%%%% Input Winding location vectors in x and y coords
astr=[ '%%Height of Winding Cross Section for each winding (in h direction)... ', num2str(numwind) ' elements are required in this vector. \n a=['] ;
% loop through the windings and designate vector element for position in
% each additional winding
for i=1:numwind
astr=[astr ' ?? '];
end
% close vector and write it to file
astr=[astr ']; %%[meters] \n \n'];
fprintf(fid,astr);
bstr=[ '%%%%Width of Winding Cross Section for each winding (in bw direction) ... ', num2str(numwind) ' elements are required in this vector. \n b=['] ;
% loop through the windings and designate vector element for position in
% each additional winding
for i=1:numwind
bstr=[bstr ' ?? '];
end
% close vector and write it to file
bstr=[bstr ']; %%[meters] \n \n'];
fprintf(fid,bstr);
%%%%% Input Winding Starting vectors in x and y coords
xzerostr=[ '%% Location of Center of Winding Cross Section for each winding (in h direction) ... ' num2str(numwind) ' elements are required in this vector. \n xzero=['] ;
% loop through the windings and designate vector element for position in
% each additional winding
for i=1:numwind
xzerostr=[xzerostr ' ?? '];
end
% close vector and write it to file
xzerostr=[xzerostr ']; %%[meters] \n \n'];
fprintf(fid,xzerostr);
yzerostr=[ '%% Location of Center of Winding Cross Section for each winding (in bw direction) ... ' num2str(numwind) ' elements are required in this vector. \n yzero=['] ;
% loop through the windings and designate vector element for position in
% each additional winding
for i=1:numwind
yzerostr=[yzerostr ' ?? '];
end
% close vector and write it to file
yzerostr=[yzerostr ']; %%[meters] \n \n'];
fprintf(fid,yzerostr);
%%%% I (current)
% loop through the windings and create matrix row to designate each additional winding
for i=1:numwind
string='%%Current at the end of the last interval is assumed equal to the first current value. \n \n';
string=[string ' %%Current values for winding ' num2str(numwind-((numwind-i))) ' at the beginning of each time segment... ' num2str(numtime) ' elements are required in this vector. \n I(' num2str(numwind-((numwind-i))) ',:)= ['];
% loop through time segments and designate vector element for each additional
% segment
for k=1:numtime
string=[string ' ?? '];
end
% close vector and write it to file
string=[string ']; %%[amperes] \n \n' ];
fprintf(fid,string);
end
%%%% N (turns)
Nstr=[ '%%Number of turns in each winding... ' num2str(numwind) ' elements are required in this vector. \n N=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
Nstr=[Nstr ' ?? '];
end
% close vector and write it to file
Nstr=[Nstr ']; %%[turns] \n \n'];
fprintf(fid,Nstr);
%%%% Plot orignal design - create vector of wire diameters
origdiastr =[ '%%Wire gauge for each winding... ' num2str(numwind) ' elements are required in this vector. \n orig_wiresize=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
origdiastr=[origdiastr ' 0 '];
end
% close vector and write it to file
origdiastr=[origdiastr ']; %%[AWG gauge] \n \n'];
fprintf(fid,origdiastr);
%%%% Plot orignal design - create vector of number of strands in each winding
orignumstrandstr =[ '%%Number of strands in each winding... ' num2str(numwind) ' elements are required in this vector. \n orig_numstrands=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
orignumstrandstr=[orignumstrandstr ' 0 '];
end
% close vector and write it to file
orignumstrandstr=[orignumstrandstr ']; %%[number of strands] \n \n'];
fprintf(fid,orignumstrandstr);
% make matlab spit out the gap/core dimensions in the file just as it would display it ordinarily
insulbuildstr=[ '\n %%The wire insulation build is \n insulbuild= ''' insulbuild '''; \n'];
fprintf(fid,insulbuildstr); % print insulation build
gaplocstr=[ '\n %%Location of Gap \n gaploc= ''' gaploc ''' ; \n'];
% gaplocstr=[ '%%Breadth of Winding Window [meters] \n bw=' num2str(bw) '; \n \n %%Height of Winding Window \n h=' num2str(h) '; \n
%%Breadth of Bobbin Window [meters] \n bb=' num2str(bb) '; \n \n %%Height of Bobbin Window \n hb=' num2str(hb) '; \n \n
%%Gaploc \n gaploc=''' gaploc '''; \n' gaplenstr];
fprintf(fid,gaplocstr); % print gap location
%this gives the user tha ablitity to edit datafile immediately
fclose(fid);
msg1=['You have created file "' (filename) '" in path "' (pathname) '" Please open this file and enter your data, then rerun litzopt and select "Existing" to calculate the optimal designs for your transformer.' ];
answer=questdlg(msg1,'Choose','Edit Now', 'Ok, I will edit it later', 'default');
if answer(1)=='E'
edit(filename);
else
cd(origdirect);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end twodsimulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% onedsimulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function is run when a new data file is to be created and the
% user selects a 1 dimensional simulation. For the analytic sol'n, the input parameters no longer
% include winding volumes and simulated field values as they are calculated
% in another function
function onedimsimulation()
prompt={'Number of windings:','Number of time segments for piecewise-linear current waveforms:'};
answer=inputdlg(prompt);
vect=str2double(answer);
numtime=vect(2);
numwind=vect(1);
uo=4e-7*pi; %#ok<*NASGU>
answer1=questdlg('Are there any gaps in the core? ','Choose','No legs gapped','Center leg gapped','All legs gapped','No legs gapped');
if answer1(1)=='C'
gaploc='Center';
elseif answer(1) == 'A'
gaploc = 'All';
else
gaploc='No Gap';
end
answer2=questdlg('What build is the wire insulation?','Choose','Single','Heavy','Single');
if answer2(1)=='S' %this is single build insulation
insulbuild='s';
else %this is heavy insulation
insulbuild = 'h';
end
% check to see if '.m' follows filename
[filename,pathname]=uiputfile('Untitled', 'Save new file as:');
both= and (filename(length(filename)-1)=='.',filename(length(filename))=='m');
% if '.m' does not follow filename, place it there
if both==0; % if false statements
filename=[filename '.m'];
end
% store original pathname and set path to new pathname
origdirect=cd;
cd(pathname);
% create new file in the desired path
[fid,errormsg]=fopen(filename,'w');
% begin writing file
startstr='%%For an example and help, go to www.thayer.dartmouth.edu/inductor/litzopt, or refer to winding data for use with litzopt. \n \n';
startstr=[startstr '%%Please replace the question marks with the pertinent data and leave all else alone....\n \n'];
startstr=[startstr '%%There are ' num2str(numwind) ' windings in your transformer. \n \n'];
startstr=[startstr '%%Resistivity (1.77e-8 is for room temperature copper) \n '];
startstr=[startstr ' rhoc = 1.77e-8; %% [ohm-meters]\n \n'];
startstr=[startstr '%%Perform optimization on designs with wire sizes ranging from American Wire Gauge strand size (32 to 48 is the default). \n'];
startstr=[startstr ' awg = [32:2:48]; %% \n \n'];
startstr=[startstr '%%Breadth of the Core Window in meters \n'];
startstr=[startstr ' bw = ?? ; %% \n \n'];
startstr=[startstr '%%Area of the Bobbin Window in meters \n'];
startstr=[startstr ' areabb = ?? ; %% \n \n'];
startstr=[startstr '%%Maximum Achieveable Packing Factor \n'];
startstr=[startstr ' userfp = ?? ; %% \n \n'];
fprintf(fid,startstr);
%%% Length, create vector of winding lengths
istr=['%%Average turn length for each winding... ' num2str(numwind) ' elements are required in this vector. \n'] ;
istr=[istr ' len =[ '];
%%%% loop through windings and designate vector element for each winding length
for i=1:numwind
istr=[istr ' ?? '];
end
% close vector and write it to file
istr=[istr ']; %%[meters] \n \n'];
fprintf(fid,istr);
%%% Time Duration %%%
% create vector of time durations for each segment
dtstr=['%%Duration of each time segment... ' num2str(numtime) ' elements are required in this vector. \n'];
dtstr=[dtstr ' dt =[ '];
% loop through number of time segments and designate a vector element for additional
% time segment duration
for j=1:numtime
dtstr=[dtstr ' ?? '];
end
% close vector and write it to file
dtstr=[dtstr ']; %%[seconds] \n \n'];
fprintf(fid,dtstr);
%%%% I (current)
% loop through the windings and create matrix row to designate each additional winding
for i=1:numwind
string='%%Current at the end of the last interval is assumed equal to the first current value. \n \n';
string=[string ' %%Current values for winding ' num2str(numwind-((numwind-i))) ' at the beginning of each time segment... ' num2str(numtime) ' elements are required in this vector. \n I(' num2str(numwind-((numwind-i))) ',:)= ['];
% loop through time segments and designate vector element for each additional
% segment
for k=1:numtime
string=[string ' ?? '];
end
% close vector and write it to file
string=[string ']; %%[amperes] \n \n' ];
fprintf(fid,string);
end
%%%% N (turns)
Nstr=[ '%%Number of turns in each winding... ' num2str(numwind) ' elements are required in this vector. \n N=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
Nstr=[Nstr ' ?? '];
end
% close vector and write it to file
Nstr=[Nstr ']; %%[turns] \n \n'];
fprintf(fid,Nstr);
%%%% Plot orignal design - create vector of wire diameters
origdiastr =[ '%%Wire gauge for each winding... ' num2str(numwind) ' elements are required in this vector. \n orig_wiresize=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
origdiastr=[origdiastr ' 0 '];
end
% close vector and write it to file
origdiastr=[origdiastr ']; %%[AWG gauge] \n \n'];
fprintf(fid,origdiastr);
%%%% Plot orignal design - create vector of number of strands in each winding
orignumstrandstr =[ '%%Number of strands in each winding... ' num2str(numwind) ' elements are required in this vector. \n orig_numstrands=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
orignumstrandstr=[orignumstrandstr ' 0 '];
end
% close vector and write it to file
orignumstrandstr=[orignumstrandstr ']; %%[number of strands] \n \n'];
fprintf(fid,orignumstrandstr);
gaplocstr=[ '%%Gap \n gap=''' gaploc '''; \n'];
fprintf(fid,gaplocstr); % print gaplocation
insulationstr=[ '%%The wire insulation build is \n insulbuild=''' insulbuild '''; \n'];
fprintf(fid,insulationstr); % print wire insulation build
fclose(fid);
msg1=['You have created file "' (filename) '" in path "' (pathname) '" Please open this file and enter your data, then rerun litzopt and select "Existing" to calculate the optimal designs for your transformer.' ];
answer=questdlg(msg1,'Choose','Edit Now', 'Ok, I will edit it later', 'default');
if answer(1)=='E'
edit(filename);
else
cd(origdirect);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%end onedimsimulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% start FEA%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function is run whenever a new data file is to be created and the
% user specifies numerical field values are to be used. The output is a
% data file with the valid matrix and vector dimensions and blanks
% (or ??'s) for the user to fill in.
function fea()
prompt={'Number of windings:','Number of time segments for piecewise-linear current waveforms:'};
answer=inputdlg(prompt);
vect=str2double(answer);
numtime=vect(2);
numwind=vect(1);
answer1=questdlg('Are there any gaps in the core? ','Choose','No legs gapped','Center leg gapped','All legs gapped','No legs gapped');
if answer1(1)=='C'
gaploc='Center';
elseif answer(1) == 'A'
gaploc = 'All';
else
gaploc='No Gap';
gaplen = 0;
end
answer2=questdlg('What build is the wire insulation?','Choose','Single','Heavy','Single');
if answer2(1)=='S' %this is single build insulation
insulbuild='s';
else %this is heavy insulation
insulbuild = 'h';
end
% check to see if '.m' follows filename
[filename,pathname]=uiputfile('Untitled', 'Save new file as:');
both= and (filename(length(filename)-1)=='.',filename(length(filename))=='m');
% if '.m' does not follow filename, place it there
if both==0; % if false statements
filename=[filename '.m'];
end
% store original pathname and set path to new pathname
origdirect=cd;
cd(pathname);
% create new file in the desired path
[fid,errormsg]=fopen(filename,'w');
% begin writing file
startstr='%%For an example and help, go to www.thayer.dartmouth.edu/inductor/litzopt, or refer to winding data for use with litzopt. \n \n';
startstr=[startstr '%%Please replace the question marks with the pertinent data and leave all else alone....\n \n'];
startstr=[startstr '%%There are ' num2str(numwind) ' windings in your transformer. \n \n'];
startstr=[startstr '%%Resistivity (1.77e-8 is for room temperature copper) \n '];
startstr=[startstr ' rhoc = 1.77e-8; %% [ohm-meters]\n \n'];
startstr=[startstr '%%Perform optimization on designs with wire sizes ranging from American Wire Gauge strand size (32 to 48 is the default). \n'];
startstr=[startstr ' awg = [32:2:48]; %% \n \n'];
startstr=[startstr '%%Breadth of the Core Window in meters \n'];
startstr=[startstr ' bw = ?? ; %% \n \n'];
startstr=[startstr '%%Height of the Core Window in meters \n'];
startstr=[startstr ' h = ?? ; %% \n \n'];
startstr=[startstr '%%Breadth of the Bobbin Window in meters \n'];
startstr=[startstr ' bb = ?? ; %% \n \n'];
startstr=[startstr '%%Height of the Bobbin Window in meters \n'];
startstr=[startstr ' hb = ?? ; %% \n \n'];
startstr=[startstr '%%Maximum Achieveable Packing Factor \n'];
startstr=[startstr ' userfp = ?? ; %% \n \n'];
fprintf(fid,startstr);
%%% Length, create vector of winding lengths
istr=['%%Average turn length for each winding...' num2str(numwind) ' elements are required in this vector. \n'] ;
istr=[istr ' len =[ '];
%%%% loop through windings and designate vector element for each winding length
for i=1:numwind
istr=[istr ' ?? '];
end
% close vector and write it to file
istr=[istr ']; %%[meters] \n \n'];
fprintf(fid,istr);
%%% Time Duration %%%
% create vector of time durations for each segment
dtstr=['%%Duration of each time segment... ' num2str(numtime) ' elements are required in this vector. \n'];
dtstr=[dtstr ' dt =[ '];
% loop through number of time segments and designate a vector element for additional
% time segment duration
for j=1:numtime
dtstr=[dtstr ' ?? '];
end
% close vector and write it to file
dtstr=[dtstr ']; %%[seconds] \n \n'];
fprintf(fid,dtstr);
%%%% I (current)
% loop through the windings and create matrix row to designate each additional winding
for i=1:numwind
string='%%Current at the end of the last interval is assumed equal to the first current value. \n \n';
string=[string ' %%Current values for winding ' num2str(numwind-((numwind-i))) ' at the beginning of each time segment... ' num2str(numtime) ' elements are required in this vector \n I(' num2str(numwind-((numwind-i))) ',:)= ['];
% loop through time segments and designate vector element for each additional
% segment
for k=1:numtime
string=[string ' ?? '];
end
% close vector and write it to file
string=[string ']; %%[amperes] \n \n' ];
fprintf(fid,string);
end
%%%% Volume
% create vector of winding volumes
volstr=['%%Volume of each winding... ' num2str(numwind) ' elements are required in this vector. \n vol=['] ;
% loop through the windings and designate vector element for volume of
% each additional winding
for i=1:numwind
volstr=[volstr ' ?? '];
end
volstr=[volstr ']; %%[cubic meters] \n \n'];
volstr=[volstr '%%Number of turns in each winding... ' num2str(numwind) ' elements are required in this vector. \n N=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
volstr=[volstr ' ?? '];
end
% close vector and write it to file
volstr=[volstr ']; %%[turns] \n \n'];
fprintf(fid,volstr);
%%% B (field)
for i= 1:numwind
bstr=['%%Integral of Bsquared values from simulation with current in winding ' num2str(i) '\n \n'];
fprintf(fid,bstr);
for j=i:numwind;
if i~=j
innerstr= ['%%Integral of B squared values from simulation with current in winding ' num2str(i) ' and ' num2str(j) '\n \n'];
fprintf(fid,innerstr);
end
for imat=1:numwind;
bigstr=' %% Enter integral of B squared ' ;
bigstr=[bigstr ' over winding ' num2str(imat) '\n' ' int_B2(' num2str(i) ',' num2str(j) ',' num2str(imat) ')= ??; %%[T^2-m^3 ] \n'];
fprintf(fid,bigstr);
end
fprintf(fid,'\n');
end
end
%%%% Plot orignal design - create vector of wire diameters
origdiastr =[ '%%Wire gauge for each winding... ' num2str(numwind) ' elements are required in this vector. \n orig_wiresize=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
origdiastr=[origdiastr ' 0 '];
end
% close vector and write it to file
origdiastr=[origdiastr ']; %%[AWG gauge] \n \n'];
fprintf(fid,origdiastr);
%%%% Plot orignal design - create vector of number of strands in each winding
orignumstrandstr =[ '%%Number of strands in each winding... ' num2str(numwind) ' elements are required in this vector. \n orig_numstrands=['] ;
% loop through the windings and designate vector element for turns in
% each additional winding
for i=1:numwind
orignumstrandstr=[orignumstrandstr ' 0 '];
end
% close vector and write it to file
orignumstrandstr=[orignumstrandstr ']; %%[number of strands] \n \n'];
fprintf(fid,orignumstrandstr);
%%%%% Input Winding location vectors in x and y coords
astr=[ '%%Height of Winding Cross Section for each winding (in h direction)... ', num2str(numwind) ' elements are required in this vector. \n a=['] ;
% loop through the windings and designate vector element for position in
% each additional winding
for i=1:numwind
astr=[astr ' ?? '];
end
% close vector and write it to file
astr=[astr ']; %%[meters] \n \n'];
fprintf(fid,astr);
bstr=[ '%%%%Width of Winding Cross Section for each winding (in bw direction) ... ', num2str(numwind) ' elements are required in this vector. \n b=['] ;
% loop through the windings and designate vector element for position in
% each additional winding
for i=1:numwind
bstr=[bstr ' ?? '];
end
% close vector and write it to file
bstr=[bstr ']; %%[meters] \n \n'];
fprintf(fid,bstr);
%%%%% Input Winding Starting vectors in x and y coords
xzerostr=[ '%% Location of Center of Winding Cross Section for each winding (in h direction) ... ' num2str(numwind) ' elements are required in this vector. \n xzero=['] ;
% loop through the windings and designate vector element for position in
% each additional winding
for i=1:numwind
xzerostr=[xzerostr ' ?? '];
end
% close vector and write it to file
xzerostr=[xzerostr ']; %%[meters] \n \n'];
fprintf(fid,xzerostr);
yzerostr=[ '%% Location of Center of Winding Cross Section for each winding (in bw direction) ... ' num2str(numwind) ' elements are required in this vector. \n yzero=['] ;
% loop through the windings and designate vector element for position in
% each additional winding
for i=1:numwind
yzerostr=[yzerostr ' ?? '];
end
% close vector and write it to file
yzerostr=[yzerostr ']; %%[meters] \n \n'];
fprintf(fid,yzerostr);
% make matlab spit out the gap location and insulation build in the file just as it would display it ordinarily
insulbuildstr=[ '\n %%The wire insulation build is \n insulbuild= ''' insulbuild '''; \n'];
fprintf(fid,insulbuildstr); % print insulation build
gaplocstr=[ '\n %%Location of Gap \n gaploc= ''' gaploc ''' ; \n'];
% gaplocstr=[ '%%Breadth of Winding Window [meters] \n bw=' num2str(bw) '; \n \n
%%Height of Winding Window \n h=' num2str(h) '; \n %%Breadth of Bobbin Window [meters] \n bb=' num2str(bb) '; \n \n
%%Height of Bobbin Window \n hb=' num2str(hb) '; \n \n%%Gaploc \n gaploc=''' gaploc '''; \n' gaplenstr];
fprintf(fid,gaplocstr); % print gap location
fclose(fid);
msg1=['You have created file "' (filename) '" in path "' (pathname) '" Please open this file and enter your data, then rerun litzopt and select "Existing" to calculate the optimal designs for your transformer.' ];
answer=questdlg(msg1,'Choose','Edit Now', 'Ok, I will edit it later', 'default');
if answer(1)=='E'
edit(filename);
else
cd(origdirect);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%end FEA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% CALCULATIONS SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% start of Existing%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function imports the data file edited by the user and decides if user-simulated
% field data has been given or the expressions derived from an analytic sol'n are to be used
function existing()
% N by M means N windings and M time segments
% user specifies file to be read
[filename, pathname] = uigetfile('Untitled');
% store original pathname and set path to new pathname
origdirect=cd;
cd(pathname);
% remove '.m' from filename and execute the file as a string
len=length(filename);
filename=filename(1:1:len-2);
eval(filename);
% Invert matrix dimensions from file so that user can input the vectors in
% horizontal fashion , N by M means N windings and M time segments
uo=4*pi*1e-7;
N=N'; %#ok<*NODEF>
len=len';
D=size(I);
orig_numstrands = orig_numstrands';
orig_wiresize = orig_wiresize';
numwind=D(1); % store the number of windings and time segments
numtime=D(2);
awgfine = min(awg)-4:.1:max(awg)+4; % Special range for figure and table
%start building new program here
I(:,numtime+1)=I(:,1); % want currents at beg of interval 1
%determine the Ims and Irms
[Ims, Irms] = Imscurrent(I, dt, numtime, numwind);
% if dimfield==2; % User has selected the 2D sol'n
% B2avgcalc2d calculates average squared field data Inputs are
% window breadth gap location, width of gap, uo, number of turns,
% current, and number of windings
% Divide by the winding volumes for FEA case prior to runnung hardcalc
if sum(strcmp(who,'int_B2')); % User has selected FEA sol'n
% set verion id to 'f' for fea sol'n
versionID = 'f';
% Divide by the winding volumes for FEA case prior to runnung avgbsquare
for pn=1:numwind
avg_B2(:,:,pn)=int_B2(:,:,pn)/vol(pn);
end
elseif sum(strcmp(who,'gaplen')); % User has selected 2D sol'n
% set verion id to 't' for 'two' 2d sol'n
versionID = 't';
%%%%%%%%% B2avgcalc2d calculates average squared field data
[avg_B2]=B2avgcalc2d(a,b,gaplen,gaploc,bw,h,xzero,yzero,xdiv,ydiv,Ximage,Yimage);
else % % User has selected 1D sol'n
% set verion id to 'o' for 'one' 1d sol'n
versionID = 'o';
%%%%%%%%% B2avgcalc1d calculates average squared field data
[avg_B2]=B2avgcalc1d(bw,gap,uo,numwind);
%calculate bb, hb,h and a from areabb, bw, N, numwind and Irms
[a,hb,bb, h] = windinght1d(Irms, bw, areabb, N, numwind);
end
% OK, avg_B2 represents area averages
% Run averagebsquare in any of the cases for the possible field geometries and get the volume averages
[B2avg]=averagebsquare(avg_B2,N);
% create current waveform figure
PWLcurrentplotL(dt, numwind, I);
%Call a function to do the calculations needed to find k_ell.
kl = find_k_ell(numwind,numtime,I,Ims,dt,N,len,rhoc,B2avg);
%Now call a function to find the actual optimal designs, plus the buildable designs (xxx2 are buildable results), plus...
%special for the paper, the full-bobbin designs (xxx3)
[cost44, totcost,totloss,loss,n,fp,totcost2,totloss2,loss2,n2,totcost3,totloss3,loss3,n3, awground, n2round,totcost2round, totloss2round, dia] = find_optimal_designs(awg,numwind,kl,Irms,len,N,rhoc,bb,a,insulbuild,userfp); %#ok<*ASGLU>
%Do it again for closely spaced wire sizes to get a smooth curve for plotting.
[totcostfine,totlossfine,totcost2fine,totloss2fine] = find_optimal_designs_curvedata(awgfine,numwind,kl,Irms,len,N,rhoc,bb,a,insulbuild,userfp);
% find cost and loss for original component design
[tot_orig_cost, tot_orig_loss] = originaldesign(orig_numstrands, orig_wiresize, Irms, len, N, rhoc, kl, cost44);
%create design labels
optdesignlabels = designlabels(dia);
bbdesignlabels = designlabels(awground);
% %Plot results in design frontiers:
plotdesignsL(tot_orig_cost, tot_orig_loss, bbdesignlabels, optdesignlabels, dia, awground, totloss, totcost, totcost2round, totloss2round, totcost3, totloss3, totcostfine,totlossfine,totcost2fine,totloss2fine)
% %plots optimal design curve with full bobbins curve, hypothetical optimum curve and best buildable curve
%plot winding window cross section
if versionID ~= 'o'
XsectionplotL(xzero, yzero, a, b, h, bw, hb);
end
%displays the optimal design table
htmltable = optdestable(optdesignlabels, totcost, totloss, numwind, n, dia, fp);
disp(htmltable);
%displays the buildable design table
bdtable = buildabletable(bbdesignlabels, totcost2round, totloss2round, numwind,n2round,awground);
disp(bdtable);
%displays original designs cost and loss if an original design was entered
if tot_orig_cost ~= 0
title = 'Original Design';
coststr = ['Relative Cost of Original Design = ' num2str(tot_orig_cost) ''];
lossstr = ['Loss of Original Design (in Watts) = ' num2str(tot_orig_loss) ''];
disp(title);
disp(coststr);
disp(lossstr);
end
%%%%%%%%%%%%%%%%%%%%%%end of existing%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%start of avergebsquare%%%%%%%%%%%%%%%%%%%%%%%%
% Generates spatial average value of the squared field with
% respect to each of the windings due to the current in all
% the windings. Inputs are 'Bsquared',the matrix containing
% square of the sum of area-averaged fields due to windings 'i' and 'j'
% (the output of a finite element analysis field approximation),
% 'vol' (winding volumes), and 'N' (turns in each winding). The output
% is a matrix of area-averaged products Bi*Bj at each location (i,j) on each page
% which is the product of fields due to windings 'i' and 'j'.
function [B2avg]=averagebsquare(Bsquared,N)
%determine the number of windings
nw=length(N);
%make a column vector of one with nw rows
one=ones(nw,1);
%make a matrix of the number of turns, i.e. first column is the number of
%turns in primary,Np, second column is the number of turns in the secondary, Ns, etc
Nn=one*N';
%Create a matrix wher Nm(1,1) = Np*Np, Nm(1,2) = Np*Ns, Nm(2,1) = Ns*Np,
Nm=Nn.*Nn';
%replace the (j,i) element with the (i,j) element
for i=1:nw
for j=1:nw
Bsquared(j,i,:)=Bsquared(i,j,:);
end
end
for pn=1:nw
%take the first page of Bsquared
B2row=Bsquared(:,:,pn);
%extracts the diagonal from B2row and makes a row vector
Bdiag=spdiags(B2row,0)';
%multiple Bdiag by one matrix to get a matrix of diagonal elements of
%B2row
Bvert=one*Bdiag;
%obtain the matrix [Bvert(1,1)+Bvert(1,1) Bvert(2,1)+Bvert(1,2); Bvert(1,2)+Bvert(2,1) Bvert(2,2)+Bvert(2,2)]
allsq.terms=Bvert+Bvert';
allsq.terms;
%first, create column vector of diagonal elements in allsq.term, then make that a square matrix with those elements on the diagonal
%and all other elements zero, then mulitple by 3/2 and a unit matrix of size nw. finally subtract this matrix from the allsq.terms matrix
sq.terms=allsq.terms-eye(nw)*diag(spdiags(allsq.terms,0),0)*(3/2);
sq.terms;
%take the matrix B2row and subtract sq.terms,, then multiple by the matrix Nm which contains the number of turns squared and divid by 2
(B2row-sq.terms); %#ok<*MNEFF>
B2avg(:,:,pn)=(B2row-sq.terms).*Nm/2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%End of avergebsquare%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%start of Costnloss%%%%%%%%%%%%%%%%%%%%%%%%%%
% Returns vectors of costs, losses, and stranding for each of the calculated
% designs where each row designates a different design.
% Inputs are 'Irms' (root mean squared
% current), 'lt' (turn length), 'N' (number of turns), 'awg' (a vector of
% strand diameters for each of the optimized designs), 'rhoc' (the resistivity
% of copper at operating temperature), and 'kl' (a parameter defined by
% time-rate of change of currents and the modified dynamic loss matrix).
% The outputs are an array called 'cost' (Cost function of litz wire is assumed
% to be proportional to '(1+(k/(DC^4))+(kk))*n*lt*N' where 'k' is 1.1e-26,
% 'kk' is 2e-9, 'n' is the number of strands in the design in question
% 'lt' is the length of a single turn, and 'N' is the number of turns.
% where, dc is strand diameter, n is strand number, lt is turn
% length, and N is turn number), 'actual_loss' (the total winding power lost
% in each design and is equal to 'Irms^2*lt*N*rhoc*Fr/n/(pi*dc^2/4)' where all
% parameters are as defined above and rhoc is the resistivity of copper at
% operating temperature), and 'n' is defined to be the number of strands implemented
% in each design.
function [cost,actual_loss,n]= costnloss(Irms,lt,kl,N,awg,rhoc)
dc = awgtomet(awg); %convert wire guage to meters for calculation
%we need to calculate the optimal values for Fr at different guages, so do that
%and put result into Fr[]
Fr=1.+1./(1.-((2*COSTm(dc))./(dc.*dCOSTm(dc)))); % gives Fr
%[awg' Fr'];
n = numstrands(dc,kl,Fr); %numstrands will solve equation (2)
%now n holds a vector of the different number of strands for optimal Fr at
%different wire guages
%now calculate actual_loss
actual_loss = (Irms^2).*lt.*N.*rhoc.*Fr./n./(pi*dc.^2/4);
% now calculate cost for the cost/loss plot
Co = 0;
Fr=1.+1./(1.-((2*COSTm(dc))./(dc.*dCOSTm(dc)))); % gives Fr
cm = COSTm(dc);
cost = (Co+(cm.*(dc.^2).*n)).*(lt*N);
%%%%%%%End of Costnloss%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%startdCostm%%%%%%%%%%%%%%%%%%%%%%%%%%%
% returns the derivative of Cm(Dc) (cost). The input is 'DC'
% which is the strand diameter
function derivative = dCOSTm(DC)
k=1.1e-26;
kk=2e-9;
derivative=-((6.*k)./(DC.^7))-((2.*kk)./(DC.^3));
%%%%%%%%%%%%%%%%End of dCostm%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%start of Numstrands%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Numstrands
% Inputs are 'DC' (strand diameter), 'kl' (a parameter independent of stranding and
% defined by the modified dynamic loss matrix and time rate of change of
% winding currents), and 'Fr' (the ratio of ac to dc resistances). The output is
% strand number for the design in question
% solves the number of strands for Fr = given values
function numberst = numstrands(DC,kl,Fr)
%Fr=1.+1./(1.-((2*COSTm(DC))./(DC.*dCOSTm(DC)))); % gives Fr
Fr1 = Fr-1 ; %our equation uses fr-1 so go ahead and do it
muo = 4e-7*pi;
[numwind,numdesigns]=size(DC);
klm = kl*ones(1,numdesigns);
numberst = sqrt(Fr1./klm./(DC.^6)*(4/pi)^3);
%%%%%%%%%End of Numstrands%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%start of awgtomet%%%%%%%%%%%%%%%%%%%%%%%
% Converts from American Wire Guage strand size measure to meters. The input is
% strand diameter as measured by the AWG standard. The output is the same
% strand diameter in units of meters.
function y = awgtomet(x)
y = 0.0050.*(92).^((36.-x)./39); %convert from awg to inches
y = (2.54e-2)*y; %convert from inches to meters
%%%%%%%%%%%%%%%%%%%%%%%%%%end of awgtomet%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% start of COSTm.m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculates cost basis proportional to the additional cost per unit mass of a
% given strand of diameter. The input is 'DC' (the strand diameter) and the
% output is cost per unit mass assuming a cost function for litz wire
% '1+(1.1e-26/(DC^6))+(2e-9/(DC^2))'.
function costunitmass = COSTm(DC)
k = 1.1e-26;
kk = 2e-9;
costunitmass = 1+(k./(DC.^6))+(kk./(DC.^2)); % calculate cost for additional litz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end of COSTm %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start of B2avgcalc1d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% B (field) calculation of area-average of B squared over each transformer winding in 1 dimension
% inputs are breadth of the winding window, gap position, permeability of free space, and number of windings
% 7/17/03 jdp switched gap = b or c in B2avgcalc1d
function [B2avg]=B2avgcalc1d(bw,gap,uo,numwind)
%the inputs:
% bw - the breadth of the winding window
% gap - the gap locations??
% uo - the permability of free space
% numwind - the number of windings
%%B=0 is on nongap side
%check to see where gap is if there is one
%%%%%% With 1 GAP ON RIGHT %%%%%%%%%%%%%%%%%%%
if gap(1)=='B'
gapnum=1; % Run the algorithm for 1 gap on right side
B1windexcbeg2avg=0; % spatial integral of squared field before enclosing any current
B1windexcwind12avg=uo^2/(3*bw^2); % spatial integral of squared field due to excited winding of one winding system
B1windexcend2avg=uo^2/bw^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B2windexcbeg2avg=0;
B2windexcwind12avg=uo^2/(3*bw^2); % spatial integral of squared field due to excited winding of one winding system
B2windexcmiddle2avg=uo^2/bw^2;
B2windexcwind22avg=7*uo^2/(3*bw^2);
B2windsexcend2avg=4*uo^2/(bw^2); % spatial integral of squared field after enclosing all current
elseif gap(1)=='C'
gapnum=2; % Run the algorithm for 1 gap on left side
B1windexcbeg2avg=uo^2/(bw^2);
B1windexcwind12avg=uo^2/(3*bw^2); % spatial integral of squared field due to excited winding of one winding system
B1windexcend2avg=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B2windexcbegavg2=4*uo^2/(bw^2); % spatial integral of squared field before enclosing any current
B2windexcwind12avg=7*uo^2/(3*bw^2); % spatial integral of squared field due to excited winding of one winding system
B2windexcmiddle2avg=uo^2/bw^2;
B2windexcwind22avg=uo^2/(3*bw^2);
B2windsexcend2avg=0; % spatial integral of squared field after enclosing all current
else %Either all or none, and we treat them the same--any excess MMF dropped equally on sides
gapnum=0; % Run the algorithm for 0 or 2 gaps
B1windexcbeg2avg=uo^2/(4*bw^2);
B1windexcwind12avg=uo^2/(12*bw^2);
B1windexcend2avg=uo^2/(4*bw^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B2windexcbeg2avg=uo^2/bw^2;
B2windexcwind12avg=uo^2/(3*bw^2); % spatial integral of squared field due to first excited (current carrying)winding in two winding system
B2windexcmiddle2avg=0;
Bmiddle2avg=0;
B2windexcwind22avg=uo^2/(3*bw^2); %spatial integral of squared field due to second excited (current carrying) winding in two windi
B2windexcend2avg=uo^2/bw^2;
end
for pn=1:numwind % define which winding we average over
for i= 1:numwind % USE 1 winding excitation
%%%% Find out where pn is wrt i
if ipn % we have not enclosed any current and have B(x)=Bbeg in our winding
B2avg(i,i,pn)=B1windexcbeg2avg;
else %i=pn and we have a ramp fcn for B(x)
B2avg(i,i,pn)=B1windexcwind12avg;
end
for j=i:numwind; % try from i+1
if i~=j % USE 2 winding excitation
%%%%%% 2 winding excitation With GAP %%%%%%%%%%%%%%%%%%%
%%% Find out where pn is wrt i and j when j is > then or = to i and i and j are distinct
allcurrent=and(pn>i,pn>j);
middle=and(pn>i,pn
else % No gap
gapmult=0;
numwindplusgap=numwind;
gaplen=0;
end
% Define gapraise
gapraise=bw/2;
% First run 'ungap' to place Ribbons of length gaplen and magnitude 'gapmult'
% and to ungap the core, then run 'imagecalc' to obtain field values
[agap,bgap,xzerogap,yzerogap]=ungap(gaplen,gapmult,gaploc,gapraise,h);
% the sequence of excitations are called called repeatedely and call 'foreachexcitation' in the pattern of winding excitations
% begin by initializing excitation vector
excwind=zeros(numwind,1);
for i= 1:numwind % USE 1 winding excitation
excwind(i)=1;
for j=i:numwind;
if i~=j % USE 2 winding excitation
excwind(i)=1;
excwind(j)=1;
end
% the fcn 'foreachexcitation' employs the method of images to generate a
% 2D field sol'n, Inputs are the full width and length of all
% windings, gap length, gap mult, breadth of core window,
% height of core window, number of turns, and coordinates of
% bottom left corner of winding cross-sections (origin at inside
% of core). The matrix of int_B2 is the output
int_B2(i,j,:)=foreachexcitation(a,b,gaplen,gapmult, xzero,yzero,bw,h,agap,bgap,xzerogap, yzerogap, xdiv,ydiv,excwind,numwindplusgap,Ximage,Yimage);
%Now the (i,j) position on all the pages has just been calculated. Store these in the paged matrix.
% And define all the (j,i) positions as well.
int_B2(j,i,:)=int_B2(i,j,:);
excwind=zeros(numwind,1);
end
end % end i loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end of B2avgcalc2d%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start of foreachexcitation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function foreachexcitation changes winding excitations in sequence to obtain all field values described
% in the SFD method. inputs are winding positions and locations in window, gap length and number, breadth and
% height of the winding window, gap locations, winding divisions for field computation, vector describing excitation and numnber
% of images to be computed
% h - height of core window
% bw - breadth of winding window
% a - Height of Winding Cross Section for each winding (in h direction)
% b - Width of Winding Cross Section for each winding (in bw direction)
% xzero - Location of Center of Winding Cross Section for each winding (in h direction)
% yzero - Location of Center of Winding Cross Section for each winding (in bw direction)
% gaplen - length of gaps
% gaploc - the location of the gaps
% xdiv - the number of divisions in the x dimension for field calc
% ydiv - the number of divisions in the y dimension for field calc
% Ximage - the number of images in the x directions for method od images
% Yimage - the number of images in the Y directions for method of images
% gapmult - a gap multiplier determine in B2avgcalc2d
% xzerogap - from function ungap
% yzerogap - from function ungap
% excwind - keeps tract of which winding is excited
% numwindplusgap - keeps tract of the number of windings and the number of
% gaps for field calc
function [int_B2]=foreachexcitation(a,b,gaplen,gapmult, xzero,yzero,bw,h,agap,bgap,xzerogap, yzerogap, xdiv,ydiv,excwind,numwindplusgap,Ximage,Yimage) %#ok<*INUSL>
disp(['Now computing field for the following pattern of excitation: ' num2str(excwind')])
numwind=length(a);
numexc = sum(excwind);
% set gap elements as excited
gaps=numwindplusgap-numwind;
for i=numwind+1:numwindplusgap %for each gap excite the ribboon with -1 current
excwind(i)=-(numexc/gaps);
end
% combine wind vect's and gap vects
a=[a agap];
b=[b bgap];
xzero=[xzero xzerogap];
yzero=[yzero yzerogap];
% this fcn computes the field due to a uniform rectangular current
% distribution, inputs are rect location, size,
[int_B2]=imagecalc(a,b,xzero,yzero,bw,h,xdiv,ydiv,excwind,numwind,Ximage,Yimage);
%%%%%%%%%%%%%%%%%%%%%%%%End for each excitation%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start of imagecalc%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function imagecalc computes all data needed to use method of images about high perbeablilty
% core, the inputs are dimension and location of each of the windings in the window, breadth and
% height of the window, number of divisions for windings, excited winding vector, and number of
% images to be computed. The entire system of imaged excited windings are generated.
function [int_B2]=imagecalc(a,b,xzero,yzero,bw,h,xdiv,ydiv,excwind,numwind,Ximage,Yimage)
% a and b are the width and height and xzero yzero are the center points
% assign excited vectors to include winding and gap locations (eliminate unexcited windings)
% define aexc,bexc, etc
aexc=a(excwind~=0);
bexc=b(excwind~=0);
xzeroexc=xzero(excwind~=0);
yzeroexc=yzero(excwind~=0);
% We will image the geometry of excited windings to obtain a ximage by yimage system and then run Rectfldcalc, each wind/rib is on a page
apage=reshape(aexc,1,1,length(aexc));
bpage=reshape(bexc,1,1,length(bexc));
xzeropage=reshape(xzeroexc,1,1,length(xzeroexc));
yzeropage=reshape(yzeroexc,1,1,length(yzeroexc));
for n=2:(Yimage-1)/2+1 % Image up for y's, leave x values alone
yabove=bw*(n-1)-yzeropage(1,:,:); % distance above winding centers to core
yzeroadd=yzeropage(1,:,:)+2.*yabove; % %define bottom row of matrix
xzeroadd=xzeropage(1,:,:);
aadd=apage(1,:,:);
badd=bpage(1,:,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yzeropage=[yzeroadd; yzeropage];
xzeropage=[xzeroadd; xzeropage];
apage=[aadd; apage];
bpage=[badd; bpage];
end
for n=2:(Yimage-1)/2+1 % Image down for y's, leave x values alone
%instead of length(yzp), use correct dimesnion--added v1.4
nyzp=size(yzeropage);
n1yzp=nyzp(1);
ybelow=bw*(n-2)+yzeropage(n1yzp,:,:); % distance above winding centers to core
yzeroadd=yzeropage(n1yzp,:,:)-2.*ybelow; % %define bottom row of matrix
xzeroadd=xzeropage(n-1,:,:);
aadd=apage(n-1,:,:);
badd=bpage(n-1,:,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yzeropage=[yzeropage; yzeroadd];
xzeropage=[xzeropage; xzeroadd];
apage=[apage; aadd];
bpage=[bpage; badd];
end
% At this point we have one column of winding groups (containing the original winding group) that is yimage high,
% it will be good to define the farthest left column that we will define as column 1 and the fathest right as ximage
% and the original as (ximage+1)/2
for n=2:(Ximage-1)/2+1
xright=h*(n-1)-xzeropage(:,n-1,:);
xzeroadd=xzeropage(:,n-1,:)+2.*xright; %
yzeroadd=yzeropage(:,n-1,:);
aadd=apage(:,n-1,:);
badd=bpage(:,n-1,:);
yzeropage=[yzeropage yzeroadd]; % % Image right for x's leave y values alone
xzeropage=[xzeropage xzeroadd];
apage=[apage aadd];
bpage=[bpage badd];
end
for n=2:(Ximage-1)/2+1
xleft=h*(n-2)+xzeropage(:,1,:);
xzeroadd=xzeropage(:,1,:)-2.*xleft; %
yzeroadd=yzeropage(:,1,:);
aadd=apage(:,1,:);
badd=bpage(:,1,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yzeropage=[yzeroadd yzeropage];
xzeropage=[xzeroadd xzeropage]; % Image left for x's, leave y values alone
apage=[aadd apage];
bpage=[badd bpage];
end
excwindonly = excwind(excwind~=0);
for pn=1:numwind
% define apn,bpn, (all single values) etc for the pnth winding in the system
apn=a(pn);
bpn=b(pn);
xzeropn=xzero(pn);
yzeropn=yzero(pn);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Compute Fields for Imaged Geometry
% run rectfldcalc over the winding pn and obtain component fields hx and hy, only pass through excited data and pn data
% hx and hy are matrices over pn that give the x and y comp of field as a fcn of pos.
[int_B2(:,:,pn)] = Rectfldcalc(apage,bpage,xzeropage,yzeropage,excwindonly,apn,bpn,xzeropn,yzeropn,bw,h,xdiv,ydiv) ;
end
%%%%%%%%%%%%%%%%%%%%%%%%End of imagecalc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start of rectfldcalc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function Rectfldcalc takes in the imaged system of excited windings and computes the area-average squared
% magnetic field over the cross-section of a single transformer winding. The inputs are the location and
% dimensions of the imaged system, the location and dimensions of the winding being averaged over, breadth and
% height of the winding window, and the number of divisions of the winding. The single area-avergae field value
% is the output.
function [avg_B2]=Rectfldcalc(aexc,bexc,xzeroexc,yzeroexc,excwind,apn,bpn,xzeropn,yzeropn,bw,h,xdiv,ydiv)
% We would like to arrive at the following for the pnth winding:
% a magnitude Hx as a function of 2 coordinates x in xspace and y in yspace in pnth wind
% a magnitude Hy as a function of 2 coordinates x in xspace and y in yspace in pnth wind
% a nd b are the width and height and xzero yzero are the center points
%keyboard
[i,j,k]=(size(aexc)); %#ok % Find number of windings/ribbons excited in system
excwindpages=ones(i,j,k);
for index=1:k
excwindpages(:,:,index)=excwind(index)*excwindpages(:,:,index);
end
halflenpn=apn./2; %define half length and width for pn winding
halfwidpn=bpn./2;
% discretize area of pn winding
ydivpn=ydiv; %=round(ydiv.*bpn/bw);
xdivpn=xdiv; %round(xdiv.*apn/h); % compute number of divisions
% find where pnth winding begins and ends
xstartpn=xzeropn-halflenpn;
ystartpn=yzeropn-halfwidpn;
yfinishpn=ystartpn+bpn;
xfinishpn=xstartpn+apn;
% Define (x,y) discretizations inside pn windings
% a and b describe center coordintes of winding/ribbons
% obtain discretization pn winding
yspace=bpn/ydivpn;
xspace=apn/xdivpn; % compute size of divisions
% n=3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n=2 ------------- (i,3)
% Obtain corner coordinates for all excited winds/ribs where xcorner(i,j) is the x position of corner j of winding i l (i,2) l
% and ycorner(i,n) is the y position of corner n of winding i l l
% l l
% l l
halfheight=aexc/2;
halfwidth=bexc/2;
halfxspace=xspace/2;
halfyspace=yspace/2;
xbotleft=xzeroexc-halfheight;
ybotleft=yzeroexc-halfwidth;
xtopleft=xzeroexc-halfheight;
ytopleft=yzeroexc+halfwidth;
xtopright=xzeroexc+halfheight;
ytopright=yzeroexc+halfwidth;
xbotright=xzeroexc+halfheight;
ybotright=yzeroexc-halfwidth;
%%%%%%%%%%%%%%%%%%%%
xstartpn=xstartpn+halfxspace; % begin at the center of the first division
ystartpn=ystartpn+halfyspace;
% Iterate through field pts by going across row, then start over on row above in pn, for every fld point get 4 r's and 4 thetas
% pn spans from xstartpn,ystartpn -------> (xfinishpn,yfinishpn)
% The values of theta are always less than 180 deg in magnitude, example 170 deg - (-170 deg) = 20 deg, not 340
twopi=2*pi;
for k=1:ydivpn
for j=1:xdivpn % go through a row of divisions inside pn, we compute radius angle and field one fld point at a time
%at the center point of the divisions
xsense = (xstartpn+(j-1)*xspace);
ysense = (ystartpn+(k-1)*yspace);
yrel = ysense - yzeroexc;
xrel = xsense - xzeroexc;
rbotleft=sqrt((xbotleft-xsense).^2+((ybotleft-ysense).^2));
rtopleft=sqrt((xtopleft-xsense).^2+((ytopleft-ysense).^2));
rtopright=sqrt((xtopright-xsense).^2+((ytopright-ysense).^2));
rbotright=sqrt((xbotright-xsense).^2+((ybotright-ysense).^2));
thetabotleft=atan2((ysense-ybotleft),(xsense-xbotleft));
thetatopleft=atan2((ysense-ytopleft),(xsense-xtopleft));
thetatopright=atan2((ysense-ytopright),(xsense-xtopright));
thetabotright=atan2((ysense-ybotright),(xsense-xbotright));
thetadifbrbl = mod(thetabotright-thetabotleft+pi,twopi)-pi;
thetadiftrtl= mod(thetatopright-thetatopleft+pi,twopi)-pi;
thetadifbltl= mod(thetabotleft-thetatopleft+pi,twopi)-pi;
thetadifbrtr=mod(thetabotright-thetatopright+pi,twopi)-pi;
% excwind contains current magnitudes of windings (which is unit excitation) and ribbons
%numimages=size(halfwidth);
%Hx=zeros(numimages(1),numimages(2),length(excwind));
%Hy=zeros(size(Hx));
Hx=excwindpages./(8*pi*halfwidth.*halfheight).*((yrel+halfwidth).*thetadifbrbl-(yrel-halfwidth).* thetadiftrtl+...
(xrel+halfheight).*log(rbotleft./rtopleft)-(xrel-halfheight).*log(rbotright./rtopright));
Hy=excwindpages./(8*pi*halfwidth.*halfheight).*((xrel+halfheight).*thetadifbltl-(xrel-halfheight).*thetadifbrtr+...
(yrel+halfwidth).*log(rbotleft./rbotright)-(yrel-halfwidth).*log(rtopleft./rtopright));
%end
Hxtot=sum(sum(sum(Hx))); % obtain x field at (j,k) position in pn due to entire system
Hytot=sum(sum(sum(Hy)));
Hxoverpn(k,j)=Hxtot;
Hyoverpn(k,j)=Hytot;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All the computations are performed for all the winding groups arranged consistently with the imaged geometry, we have computed the
% x and y fields at each of the divisions of pn due to all the windings in all the winding groups and stored the total resulting field values
% in Hxoverpn and Hyoverpn, we now have a couple of 2 dim matrix which spans pn and x and y fields at each division
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Obtain int_B2 by adding up area squared-field products
dA=xspace*yspace; % incremental area inside winding pn
u0 = 4e-7*pi;
Bx=u0*Hxoverpn;
By=u0*Hyoverpn;
Bsquared=(Bx.^2+By.^2) ;
%Bsquared=B.^2; % obtain Bsquared for each increment of area inside pn
BsquareddA=dA.*Bsquared;
int_B2=sum(sum(BsquareddA));
avg_B2=int_B2/(apn*bpn);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end rectfldcalc%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function ungap%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [agap,bgap,xzerogap,yzerogap]=ungap(gaplen,gapmult,gaploc,gapraise,h)
%this function removes gaps in a design and repaces with a ribbon of
%current of field calc. All the inputs were calculated in B2avgcalc2d where function is called.
if gapmult==1 % we have 1 center leg gap (amend a,b,xzero,yzero)
bgap=gaplen;
agap=1e-6; % width of gap is 0
yzerogap=gapraise;
xzerogap=0;
elseif gapmult==.5; % we have 2 gaps in the outer legs only consider the first for our window
bgap=[gaplen gaplen]; % width of gap is 0
yzerogap=[gapraise gapraise];
agap=[1e-6 1e-6];
xzerogap=[0 h];
else agap=[];
bgap=[];
xzerogap=[];
yzerogap=[];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%end ungap%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%started adding new functions here 6/23/03%%%%%
%this function plots the current waveform for the local version of litzopt
function f1 = PWLcurrentplotL(dt, numwind, I) %#ok<*STOUT>
%plots current waveform
%the input:
% I - the amplitude of the waveform at the times given in the dt vector
% dt - the vector of time segments
% numwind - the number of windings
t = [0 cumsum(dt)];
figure('name','Current Waveforms')
for i=1:numwind
%do 'numwind' Current vs time plots above each other
subplot(numwind,1,i), plot(t,I(i,:),'linewidth',2.5);
ylbl=['Winding ' num2str(i)];
ylabel(ylbl, 'fontsize', 12);
grid;
hold on;
end
subplot(numwind,1,numwind),xlabel('time [sec]', 'fontsize', 12);
subplot(numwind,1,1),title('Current Waveform for Each Winding [Amps]', 'fontsize', 14);
hold off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end PWLcurrentplotL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start imscurrent%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ims, Irms] = Imscurrent(I, dt, numtime, numwind)
% compute mean-square current Ims and root mean squared current Irms
%the input:
% I - the amplitude of the waveform at the times given in the dt vector
% dt - the vector of time segments
% numtime - the number of time segments
% numwind - the number of windings
Ims = zeros(numwind,1);
for i = 1:numtime
Ims = Ims + ((I(:,i+1)-I(:,i)).^2./3 + (I(:,i)+I(:,i+1)).^2)./4./sum(dt).*dt(i);
end
Irms = sqrt(Ims) %#ok % Irms is obtained from the square root of the means square current
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end Imscurrent%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start find_k_ell%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function kl = find_k_ell(numwind,numtime,I,Ims,dt,N,len,rhoc,B2avg)
%calculate di/dt for each winding and approximate loss by SFD method
%the input:
% I - the amplitude of the waveform at the times given in the dt vector
% dt - the vector of time segments
% numtime - the number of time segments
% numwind - the number of windings
% Ims - the mean squared current
% N - the number of turns
% len - the average length of a turn
% rhoc - the resistivity of copper
% B2avg - the avgerage of the squared field
avgcont = zeros(numwind,numwind);
% calculate time rate of change of current for all time segments
% and add up contibutions and store their spatial average into 'avgcont'
%dv = (I(2)-I(1))/dt(1)
for i = 1:numtime
dv = (I(:,i+1)-I(:,i))/dt(i);
didt2 = dv*dv';
lm = didt2*dt(i)/sum(dt);
avgcont = avgcont+(lm);
end
num=zeros(size(len));
% Find Modified Dynamic Loss Matrix describing each winding
for i=1:numwind
gammahat(i)=N(i)*len(i)/pi/rhoc/4;
Dhat(:,:,i) = gammahat(i)* B2avg(:,:,i);
matrix(:,:,i) = (avgcont.*Dhat(:,:,i));
num(i)=sum(sum(matrix(:,:,i)));
end
%note: num is this number [di1/dt di2/dt] D_jhat [di1/dt; di2/dt]
% kl is defined by [di1/dt di2/dt] D_jhat [di1/dt]
% [di2/dt]
% -------------------------------
% I_rms^2 l_w rhoc
% we know that l_t*N=l_w
kl = num./(Ims.*len.*N*rhoc);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%rnd find_k_ell%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start find optimal designs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5555
function [cost44, totcost,totloss,lossW,ncorrect,fp,totcost2,totloss2,loss2,n2,totcost3,totloss3,loss3,n3, awground, n2round,totcost2round, totloss2round, dia] = find_optimal_designs(awg,numwind,kl,Irms,len,N,rhoc,bb,a,insulbuild,userfp)
%find the designs for wire sized in vector awg, and also the set of full-bobbin designs (all the "2" ones)
%numwind, kl, Irms, len, N, rhoc are all pretty obvious.
%userfp is the maximum packing factor relative to perfect square packing of cylinders that the user specifies; used for finding
%full-bobbin designs.
%the input:
% awg - a vector of wire sizes
% kl - the value of determined in function find_k_ell
% usefp - the user specificed packing factor
% numwind - the number of windings
% Irms - the root mean squared current
% N - the number of turns
% len - the average length of a turn
% rhoc - the resistivity of copper
% insulbuild - specifices if the insulation build on the wire is single or heavy build
% bb - breadth of bobbin window in meters
% a - a vector that conatins the height of each winding in meters
%CRS, 11/15/02
%JDP implement totaldesign function 1/31/03
%JDP implemented correct version of determining optimal designs with different size wire for each winding
%function called by this function
%Findtable, fprime, fprimeinverse, ff, dlookup, stranding, optpactfactor, buildabledesign, totaldesigns, costnloss2
%3/25/03 added function dlookupinverse and totaldesignAWG
%7/30/03 JDP - added roundtoAWGrange function to round dia matrix to the available AWG sizes given in awg vector
numdesigns=length(awg);
costW=zeros(numwind,numdesigns);
lossW = costW;
dia = costW;
ncorrect = costW;
awg44 = 44; % define awg44 for use i determining the kcost44, kloss44
awgfine = min(awg)-4:.1:max(awg)+4; % Special range for figure
%determine winding with highest Irms*N product
[IrmsNproduct] = Irms.*N;
[maxindex, masterwinding] = max(IrmsNproduct);
%call to costnloss determines kcost and kloss which equal the cost and loss for #44 wire
for i=1:numwind
if isfinite(kl(i)) ==1
[kcost(i,:),kloss(i,:), nk(i,:)]= costnloss(Irms(i),len(i),kl(i),N(i),awg44,rhoc);
end
end
%determine mastercost, masterloss, and mastern for masterwinding with highest Irms*N product
%call to costnloss fills the look up table with awgfine
i=masterwinding;
if isfinite(kl(i)) ==1
[mastercost, masterloss, mastern]=costnloss(Irms(i),len(i),kl(i),N(i),awg,rhoc);
[table_cost, table_loss,table_n]=costnloss(Irms(i),len(i),kl(i),N(i),awgfine,rhoc);
end
% set up emtpy matices
costW = zeros(size(mastercost));
lossW = zeros(size(masterloss));
dia = awg;
%start building the ncorrect matrix with stranding for masterwinding
ncorrect = mastern;
%find the table to look up info
table=findtable(kcost,kloss,masterwinding,table_cost,table_loss,awgfine);
%
% %this code prints out the lookup table for debugging
% table.deriv = table.deriv';
% table.cost4deriv = table.cost4deriv';
% table.cost = table.cost';
% table.loss = table.loss';
% table.awg = table.awg';
% table.deriv(161,1) = table.deriv(160,1);
% table.cost4deriv(161,1) = table.cost4deriv(160,1);
% lookuptable = [table.deriv table.cost4deriv table.cost table.loss table.awg];
kcostmatrix=(kcost*ones(1,numdesigns));
klossmatrix=(kloss*ones(1,numdesigns));
firstmatrix = ( (kloss(masterwinding,:)./kloss).*(kcost/kcost(masterwinding,:)) )*fprime((mastercost/kcost(masterwinding,:)),table);
%test firstmatrix for NaN
[firstmatrix, kcostmatrix, klossmatrix] = testnan(firstmatrix, kcostmatrix, klossmatrix);
costW = (kcostmatrix).*(fprimeinverse(firstmatrix,table));
%test costW for NaN
[costW, kcostmatrix, klossmatrix] = testnan(costW, kcostmatrix, klossmatrix);
lossW = klossmatrix.*(ff(costW./kcostmatrix,table));
%test lossW for NaN
[lossW, costW, kcostmatrix, klossmatrix] = testnanloss(lossW, costW, kcostmatrix, klossmatrix);
dia = dlookup((costW./kcostmatrix),table);
%test costW for NaN
[dia, costW, lossW, kcostmatrix, klossmatrix] = testnandia(dia, costW, lossW, kcostmatrix, klossmatrix);
%determine ncorrect, the correct number of stranding for each windings and design
ncorrect = stranding(Irms,len,kl,N,rhoc,dia,lossW);
% round the dia matrix, which is the matrix of the wire gauge sizes for each winding to AWG sizes in range
[diaround] = roundtoAWGrange(dia, awg);
% find the cost for the rounded designs
[row,numdesignround] = size(diaround);
kcostmatrix = kcost*ones(1,numdesignround);
costWroundD = (kcostmatrix).*(dlookupinverse(diaround,table));
%find the losses for the rounded wire guage matrix
klossmatrix = kloss*ones(1, numdesignround);
lossWroundD = klossmatrix.*ff(costWroundD./kcostmatrix,table);
%calculate the stranding with the diaround, the matrix of rounded diameters and loss matrix
ncorrectALL = stranding(Irms,len,kl,N,rhoc,diaround,lossWroundD);
%determine cost44 for determing relative cost
cost44=sum(kcost);
%this determines the packing factor
fp = optpackfactor(diaround, N, ncorrectALL, bb, a, insulbuild);
%Now do the buildable designs
n2 = buildabledesign(diaround, a, ncorrectALL, fp, N, userfp, bb, insulbuild);
%the totaldesigns function rounds the number of strands in each winding
[awground, n2round] = totaldesigns(n2, diaround, numwind);
% % find the cost for the final designs with rounded awg and stranding
[cost2,loss2]= costnloss2(Irms,len,n2round,N,awground,rhoc,kl);
%Now do the full-bobbin designs
n3 = fullbobbindesigns(awg, a, ncorrect, fp, N, userfp, bb, insulbuild);
%loop through different windings to obtain cost and loss vectors for full bobbin designs
for i=1:numwind
if isfinite(kl(i)) ==1
[cost3(i,:),loss3(i,:)]= costnloss2(Irms(i),len(i),n3(i,:),N(i),awg,rhoc,kl(i));
end
end
if numwind>1
totcost=sum(costW)/cost44; %total cost for optimal design
totloss=sum(lossW); %total loss for optimal design
totcost2=sum(cost2)/cost44; %total cost for buildable design
totloss2=sum(loss2); %total loss for buildable design
totcost2round=sum(cost2)/cost44; %total cost for ROUNDED buildable design
totloss2round=sum(loss2); %total loss for ROUNDED buildable design
totcost3=sum(cost3)/cost44; %total cost for buildable design
totloss3=sum(loss3); %total loss for buildable design
else
totcost=costW/cost44;
totloss=lossW;
totcost2=cost2/cost44;
totloss2=loss2;
totcost2round=cost2/cost44;
totloss2round=loss2;
totcost3=cost3/cost44;
totloss3=loss3;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end find optimal designs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start find potimal design curvedata%%%%%%%%%%%%%%%%%%%%%%%%
function [totcostfine,totlossfine,totcost2fine,totloss2fine] = find_optimal_designs_curvedata(awg,numwind,kl,Irms,len,N,rhoc,bb,a,insulbuild,userfp)
%find the designs for wire sized in vector awg, and also the set of full-bobbin designs (all the "2" ones)
%numwind, kl, Irms, len, N, rhoc are all pretty obvious.
%userfp is the maximum packing factor relative to perfect square packing of cylinders that the user specifies; used for finding full-bobbin designs.
%the inputs:
% awg - a vector of wire sizes
% kl - the value of determined in function find_k_ell
% usefp - the user specificed packing factor
% numwind - the number of windings
% Irms - the root mean squared current
% N - the number of turns
% len - the average length of a turn
% rhoc - the resistivity of copper
% insulbuild - specifices if the insulation build on the wire is single or heavy build
% bb - breadth of bobbin window in meters
% a - a vector that conatins the height of each winding in meters
%JDP implement totaldesign function 1/31/03
%JDP implemented correct version of determining optimal designs with different size wire for each winding
%function called by this function
%Findtable, fprime, fprimeinverse, ff, dlookup, stranding, optpactfactor, buildabledesign, totaldesigns, costnloss2
%3/25/03 added function dlookupinverse and totaldesignAWG
numdesigns=length(awg);
costW=zeros(numwind,numdesigns);
lossW = costW;
dia = costW;
ncorrect = costW;
awg44 = 44; % define awg44 for use i determining the kcost44, kloss44
awgfine = min(awg):.1:max(awg); % Special range for figure
%determine winding with highest Irms*N product
[IrmsNproduct] = Irms.*N;
[maxindex, masterwinding] = max(IrmsNproduct);
%call to costnloss determines kcost and kloss which equal the cost/loss for #44 wire
for i=1:numwind
if isfinite(kl(i)) ==1
[kcost(i,:),kloss(i,:), nk(i,:)]= costnloss(Irms(i),len(i),kl(i),N(i),awg44,rhoc);
end
end
%determine mastercost, masterloss, and mastern for masterwinding with highest Irms*N product
%call to costnloss fills the look up table with awgfine
i=masterwinding;
if isfinite(kl(i)) ==1
[mastercost, masterloss, mastern]=costnloss(Irms(i),len(i),kl(i),N(i),awg,rhoc);
[table_cost, table_loss, table_n]=costnloss(Irms(i),len(i),kl(i),N(i),awgfine,rhoc);
end
% set up emtpy matices
costW = zeros(size(mastercost));
lossW = zeros(size(masterloss));
dia = awg;
%start building the ncorrect matrix with stranding for masterwinding
ncorrect = mastern;
%find the table to look up info
table=findtable(kcost,kloss,masterwinding,table_cost,table_loss,awgfine);
% %this code prints out the lookup table for debugging
% table.deriv = table.deriv';
% table.cost4deriv = table.cost4deriv';
% table.cost = table.cost';
% table.loss = table.loss';
% table.awg = table.awg';
% table.deriv(161,1) = table.deriv(160,1);
% table.cost4deriv(161,1) = table.cost4deriv(160,1);
% lookuptable = [table.deriv table.cost4deriv table.cost table.loss table.awg];
kcostmatrix=(kcost*ones(1,numdesigns));
klossmatrix=(kloss*ones(1,numdesigns));
%calculate the normalized cost, loss and diameter for rest of windings
firstmatrix = ( (kloss(masterwinding,:)./kloss).*(kcost/kcost(masterwinding,:)) )*fprime((mastercost/kcost(masterwinding,:)),table);
costW = (kcostmatrix).*fprimeinverse(firstmatrix,table);
lossW = klossmatrix.*ff(costW./kcostmatrix,table);
dia = dlookup(costW./kcostmatrix,table);
%determine ncorrect, the correct number of stranding for each windings and design
ncorrect = stranding(Irms,len,kl,N,rhoc,dia,lossW);
% round the dia, the matrix of the wire gauge size for each winding
[diaround] = totaldesignsAWG(dia, awg);
% find the cost for the rounded designs
[row,numdesignround] = size(dia);
kcostmatrix = kcost*ones(1,numdesignround);
costWroundD = (kcostmatrix).*dlookupinverse(dia,table);
%find the losses for the rounded wire guage matrix
klossmatrix = kloss*ones(1, numdesignround);
lossWroundD = klossmatrix.*ff(costWroundD./kcostmatrix,table);
%calculate the stranding with the diaround, the matrix of rounded diameters and loss matrix
ncorrectALL = stranding(Irms,len,kl,N,rhoc,diaround,lossWroundD);
%determine cost44 for determing relative cost
cost44=sum(kcost);
%this determines the packing factor
fp = optpackfactor(dia, N, ncorrect, bb, a, insulbuild);
%Now do the buildable designs
n2 = buildabledesign(dia, a, ncorrect, fp, N, userfp, bb, insulbuild);
% % find the cost for the final designs with rounded awg and stranding
[cost2,loss2]= costnloss2(Irms,len,n2,N,dia,rhoc,kl);
%Now do the full-bobbin designs
n3 = fullbobbindesigns(awg, a, ncorrect, fp, N, userfp, bb, insulbuild);
%loop through different windings to obtain cost and loss vectors for full bobbin designs
for i=1:numwind
if isfinite(kl(i)) ==1
[cost3(i,:),loss3(i,:)]= costnloss2(Irms(i),len(i),n3(i,:),N(i),awg,rhoc,kl(i));
end
end
if numwind>1
totcostfine=sum(costW)/cost44; %total cost for optimal design
totlossfine=sum(lossW); %total loss for optimal design
totcost2fine=sum(cost2)/cost44; %total cost for buildable design
totloss2fine=sum(loss2); %total loss for buildable design\
totcost2round=sum(cost2)/cost44; %total cost for ROUNDED buildable design
totloss2round=sum(loss2); %total loss for ROUNDED buildable design
else
totcostfine=costW/cost44;
totlossfine=lossW;
totcost2fine=cost2/cost44;
totloss2fine=loss2;
totcost2roundfine=cost2/cost44;
totloss2roundfine=loss2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end find optmal design curvedata%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start design labels%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
function designlabels = designlabels(dia)
%this function creates the labels for the designs
% built by JDP on 5/28/03
%the inputs:
% dia - a vector of diameters
[row, design] = size(dia);
count = 0;
designlabels = [];
d = {'d'};
for i = 1:design
count = count+1;
label=num2str(count);
templabel = strcat(d,label);
designlabels = [designlabels; templabel];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end design labels%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start plotdesignsL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f2 = plotdesignsL(tot_orig_cost, tot_orig_loss, bbdesignlabels, optdesignlabels, dia, awground, totloss, totcost, totcost2round, totloss2round, totcost3, totloss3, totcostfine,totlossfine,totcost2fine,totloss2fine) %#ok<*INUSD>
%Plot results:plots optimal design frontier
% tot_orig_cost - the relative cost for the original winding design
% tot_orig_loss - the total loss for the original winding design
% bbdesignlabels - the design labels for the buildable design
% optdesignlabels - the design labels for the optimal designs
% dia - a vector fo diameters
% awground - a vector of rounded awg sizes
% totloss - the total loss for hypotheitcal optimum design - rough data
% totcost - the realtive cost for hypotheitcal optimum design - rough data
% totcost2round - the relative cost for the buildable designs
% totloss2round - the total loss for the buildable designs
% totcost3 - the relative cost of the full bobbin designs
% totloss3 - the total loss for the full bobbin designs
% totcostfine - a vector of cost results for hypothetical designs with great resolution - used for the line
% totlossfine - a vector of loss results for hypothetical designs with great resolution - used for the line
% totcost2fine - a vector that plots the cost results for buildable designs with great resolution
% totloss2fine - a vector that plots the loss results for buildable designs with great resolution
[row, numdesigns] = size(dia);
[rowbb, numbbdesigns] = size(awground);
figure('name','Optimal Design Frontier')
% totcostfine = totcostfine;
% totcost = totcost;
% totcost2round = totcost2round;
% tot_orig_cost = tot_orig_cost;
%use_opt = logical(totloss == totloss2); %where the bobbin is underfilled
% %with fine data for smooth curves
% %first plot optimal frontier in red; will be color for unbuildable designs
loglog(totcostfine,totlossfine,'-.r','linewidth',1.5)
hold on
% %Now coarse data with markers
loglog(totcost,totloss,'r^') %plots over-full bobbin designs
%plots orignal design on optimal design frontier
if tot_orig_cost ~= 0
loglog(tot_orig_cost, tot_orig_loss,'b*')
end
%if you want to plot the full bossin designs, use this code
% % %now full-bobbin designs in blue
% % %%%%%%%%%%%%%%%%%%%
% loglog(totcost3,totloss3,'--b','linewidth',1.5)
% % %%%%%%%%%%%%%%%%%%%%%
% loglog(totcost3,totloss3,'bs') %plots full bobbin designs
%and finally, the buildable designs in black
%loglog(totcost2fine,totloss2fine,'-k','linewidth',1.5)
loglog(totcost2round,totloss2round,'-k','linewidth',1.5)
%Now coarse data for markers
loglog(totcost2round,totloss2round,'ko') %plots buildable design
%label the designs on the plot
for j= 1:numbbdesigns
text(totcost2round(j)*1.08,totloss2round(j)*1.08,bbdesignlabels(j))
%text(totcost3(j)*1.05,totloss3(j)*1.05,num2str(dia(1,j)'),'color','b')
end
ylabel('Loss (W)', 'fontsize', 12),xlabel('Relative Cost', 'fontsize', 12)
set(gca,'fontsize',12)
if tot_orig_cost == 0
legend('Hypothetical optimal designs','','Best buildable designs','');
else
legend('Hypothetical optimal designs','Original design','Best buildable designs','');
end
grid
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end plotdesignsL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start Xsectionplot%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f3 = XsectionplotL(xzero, yzero, a, b, h, bw,hb)
% updated by JDP 6/24/03 to fill rectangles and label better
% this function returns a plot of the winding window
% here are the definitions of h,bw,hb,a,b,xzero and yzero for reference
% h - height of core window
% bw - breadth of winding window
% hb - height of bobbin window
% a - Height of Winding Cross Section for each winding (in h direction)
% b - Width of Winding Cross Section for each winding (in bw direction)
% xzero - Location of Center of Winding Cross Section for each winding (in h direction)
% yzero - Location of Center of Winding Cross Section for each winding (in bw direction)
figure('name','Winding Window Cross Section');
xzero = xzero.*1000; %converts from m to mm
yzero= yzero.*1000; %converts from m to mm
a = a.*1000; %converts from m to mm
b = b.*1000; %converts from m to mm
hold off
x(1,:)=[1 0 0]; %red- first winding
x(2,:)=[0 1 0]; %green- second winding
x(3,:)=[0 0 1]; %blue- third winding
x(4,:)=[1 1 0]; %yellow- fourth winding
x(5,:)=[1 0 1]; %purplish- fifth winding
xtotal = zeros(length(a),3);
%create color vector of length numwind
for g = 1:length(a)
%find the modulus of index wrt 5
mod5 = mod(g,5);
%test mod5
if mod5 == 0
xtotal(g,:) = x(5,:);
elseif mod5 == 4
xtotal(g,:) = x(4,:);
elseif mod5 == 3
xtotal(g,:) = x(3,:);
elseif mod5 == 2
xtotal(g,:) = x(2,:);
else
xtotal(g,:) = x(1,:);
end
end
%plot windings cross section and fill rectangle
for winding=1:length(a)
rectangle('position',[(yzero(winding)-b(winding)/2),(xzero(winding)-a(winding)/2),b(winding),a(winding)], 'EdgeColor', xtotal(winding,:))
fillx = [(yzero(winding)-b(winding)/2),(yzero(winding)-b(winding)/2+b(winding)),(yzero(winding)-b(winding)/2+b(winding)),(yzero(winding)-b(winding)/2)];
filly =[(xzero(winding)-a(winding)/2),(xzero(winding)-a(winding)/2),(xzero(winding)-a(winding)/2+a(winding)),(xzero(winding)-a(winding)/2+a(winding))];
fill(fillx, filly, xtotal(winding,:))
hold on
end
rectangle('position',[0,0,(bw.*1000),(hb.*1000)],'linewidth',2)
hold off
hold on
axis equal
xlabel('bw-- Window Breadth (mm)', 'fontsize', 14)
ylabel('h-- Window Height (mm)', 'fontsize', 14)
%create legend
for j=1:length(a)
title = ['Winding ' num2str(j) ' '];
text(yzero(j),xzero(j),title,'horizontalalignment','center')
end
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end XsectionplotL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start findtable%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function table=findtable(kcost,kloss,masterwinding,table_cost,table_loss,awgfine)
%this function find the lookup table for determining the optimal design where each winding
%has different size wire gauge. It normalizes cost and loss, then find the derivative of the cost function
% kcost - normalized cost constant
% kloss - normalized loss constant
% masterwinding - tells which winding is the mater winding
% table_cost - the cost for the table - the same size as awgfine
% table_loss - the loss for the lookup table - the same size as awgfine
% awgfine - a vector of wire sizes in 0.1 increments
%find normqlized cost $ loss
normalcost = table_cost./kcost(masterwinding);
normalloss = table_loss./kloss(masterwinding);
%find derivative of cost function
[deriv, cost4deriv] = numderiv(normalcost, normalloss);
table.deriv = abs(deriv);
table.cost4deriv=cost4deriv;
table.cost=normalcost;
table.loss=normalloss;
table.awg = awgfine;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end findtable%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555start numderiv%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [deriv, location] = numderiv(x,y)
% this function finds the numerical derivative
% x and y can be vectors
% created 2/28/03 by Jenna Pollock
deltax = diff(x);
deltay = diff(y);
deriv = abs(deltay)./deltax;
location= (x(1:(length(x)-1)) + x(2:length(x))/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end numderiv%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start fprime%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function deriv=fprime(normcost,table)
%this function looks up the derivative in the llokup table at the the value of normcost
%the inputs:
%normcost - the normalized cost
%table - the lookup table
deriv = interp1(table.cost4deriv,table.deriv,normcost);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end fprime%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555start fprimeinverse%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function normalcost = fprimeinverse(deriv,table)
%this function finds the normalized cost of given a derivative and the lookup table
%the inputs:
%deriv - the derviative under consideration
%table - the lookup table
normalcost = interp1(table.deriv,table.cost4deriv,deriv);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end fprimeinverse%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start dlookup%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function awg = dlookup(normcost,table)
%this function looks up the diameter cooresponding to a normalized cost in the look up table.
%the inputs:
%normcost - the normalized cost
%table - the lookup table
awg = interp1(table.cost,table.awg,normcost);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end dlookup%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start dlookupinverse%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function normcost = dlookupinverse(diaround,table)
%this function looks up the normalized cost cooresponding to a diameter(awg) in the look up table.
%the inputs:
%diaround - the rounded diameters
%table - the lookup table
normcost = interp1(table.awg,table.cost,diaround);
%%%%%%%%%%%%%%%%%%%%%%%%%%%end dlookupinverse%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start ff%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function normloss=ff(normcost,table)
%this function is the forward function
%given a noralized cost, it looksup the normalized loss, normloss
%the inputs:
%normcost - the normalized cost
%table - the lookup table
normloss = interp1(table.cost,table.loss,normcost);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end ff%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start stranding%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stranding function
% Returns vectors of costs, losses, and stranding for each of the calculated
% designs where each row designates a different design.
% Inputs are
% 'Irms' (root mean squared current),
% 'lt' (turn length),
% 'N' (number of turns),
% 'awg' (a vector of strand diameters for each of the optimized designs),
% 'rhoc' (the resistivity of copper at operating temperature)
% 'kl' (a parameter defined by time-rate of change of currents and the modified dynamic loss matrix)
%
% The outputs:
% 'n' is defined to be the number of strands implemented in each design.
function [n]= stranding(Irms,lt,kl,N,rhoc,d2,lossW2)
dc = awgtomet(d2); %convert wire guage to meters for calculation
%we need to calculate the optimal values for Fr at different guages, so do that
%and put result into Fr[]
Fr=1.+1./(1.-((2*COSTm(dc))./(dc.*dCOSTm(dc)))); % gives Fr
n = numstrands(dc,kl,Fr); %numstrands will solve equation (2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%End of stranding%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start of totaldesignsAWG%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [diaround] = totaldesignsAWG(dia, awg)
% this function will determine the diameter for all the rounded designs
% the criterion are for numbers with a remainer < 0.2, it is rounded down.
% If the remainder falls between 0.2 and 0.8, it is rounded both ways
% if the remainder is >0.8 the number id rounded up
%the inputs:
%awg - a vector of wire sizes
%diameters - vector of wire diameters
[numwind, numdesigns] = size(dia);
diamatrix = zeros(numwind,1);
tolerance = 0.4;
max_to_round = 100;
diafinal = [];
diaout =[];
awgout=[];
numbernew = 1;
%loop thru each wire size
for i = 1:numdesigns %wire size
numbernew=0;
dialow = floor(dia(:,i));
remainder = dia(:,i) - dialow;
remain = remainder;
%prune remain
in_middle = ((remainder >= tolerance)&(remainder <= (1-tolerance)));
low_enough_to_care = dia(:,i) < max_to_round;
high_enough = dia(:,i) >1;
to_keep = in_middle & low_enough_to_care & high_enough;
remain = remain(to_keep);
diasize = size(remain);
%end of pruning
%round all down, or something like that
diamatrix = dialow + ((remainder > (1-tolerance))|(dialow<1));
diafinal = [diafinal diamatrix];
diasize = size(diafinal);
numbernew = numbernew+1;
%Done round all down step
for k = 1:( length(remain))
[maxremain, index] = max(remain);
diamatrix = dialow + ((remainder >= maxremain)|(dialow<1));
diafinal = [diafinal diamatrix];
diasize = size(diafinal);
remain(index) = 0;
numbernew = numbernew + 1;
end
awgout = [awgout, awg(i)*ones(1,numbernew)];
end
diaround= diafinal;
awground = awgout;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end totaldesignAWG%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start optpackfactor%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fp] = optpackfactor(awg, N, n, bb, a, insulbuild)
%this function will determine the packing factor,fp of optimal the optimal design considering
%the insulation build
%the inputs:
% awg - a vector of awg wires sizes
% N - a vector of the number of turns in each winding
% n - a vector of the number of strands in each design
% bb - the breadth of the bobbin window
% a - a vector containing the height of each winding
% insulbuild - the insulation build selected by user, either single or
% heavy
[numwind, numdesigns]=size(awg);
fp=zeros(numwind,numdesigns);
if insulbuild == 's'
diameter = singlebuild(awg);
else
diameter = heavybuild(awg);
end
for i=1:numwind
windingarea = bb*a(i);
for j=1:numdesigns
fp(i,j) = (diameter(i,j)^2)*n(i,j)*N(i)/windingarea;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5end optpackfactor%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start buildabledesigns%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [n2] = buildabledesign(awg, a, n, fp, N, userfp, bb, insulbuild)
%this function will determine the a buildable design considering the
%insulation build. It returns the number of strands, n2, that can fit in a
%buildable design for each design
%the inputs:
% awg - a vector of awg wires sizes
% n - a vector of the number of strands in each design
% fp - the packing factor required to build an optimal design
% N - a vector of the number of turns in each winding
% userfp - the packing factor specified by the user
% bb - the breadth of the bobbin window
% a - a vector containing the height of each winding
% insulbuild - the insulation build selected by user, either single or
% heavy
[numwind, numdesigns]=size(awg);
if insulbuild == 's'
diameter = singlebuild(awg);
else
diameter = heavybuild(awg);
end
for i=1:numwind
for j=1:numdesigns
if fp(i,j)>userfp
n2(i,j) = userfp*bb*a(i)/((diameter(i,j)^2)*N(i));
else
n2(i,j) = n(i,j);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end buildable designs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start totaldesigns%%%%%%%%%%%%%%%%%%%%%%%%5
function [awground, n2round] = totaldesigns(n2, awg, numwind)
% this function will determine all the rounded designs
% the criterion are for numbers with a remainer < 0.2, it is rounded down.
% If the remainder falls between 0.2 and 0.8, it is rounded both ways
% if the remainder is >0.8 the number id rounded up
%the inputs
% n2 - a vector of the number of strands in each buildable design
% awg - a vector of wire awg sizes
% numwind - the number of windings in the component
[numwind, numdesigns] = size(awg);
n = zeros(numwind,1);
tolerance = 0.4;
max_to_round = 10;
nfinal = [];
nout =[];
awgout=[];
numbernew = 1;
%loop thru each wire size
for i = 1:numdesigns %wire size
numbernew=0;
nlow = floor(n2(:,i));
%nout = nlow;
remainder = n2(:,i) - nlow;
remain = remainder;
%prune remain
in_middle = ((remainder >= tolerance)&(remainder <= (1-tolerance)));
low_enough_to_care = n2(:,i) < max_to_round;
high_enough = n2(:,i) >1;
to_keep = in_middle & low_enough_to_care & high_enough;
remain = remain(to_keep);
nsize = size(remain);
%end of pruning
%round all down, or something like that
n = nlow + ((remainder > (1-tolerance))|(nlow<1));
nfinal = [nfinal n];
nsize = size(nfinal);
numbernew = numbernew+1;
%Done round all down step
for k = 1:( length(remain))
[maxremain, index] = max(remain);
n = nlow + ((remainder >= maxremain)|(nlow<1));
nfinal = [nfinal n];
nsize = size(nfinal);
remain(index) = 0;
numbernew = numbernew + 1;
end
awgout = [awgout, awg(:,i)*ones(1,numbernew)];
end
n2round= nfinal;
awground = awgout;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end totaldesign%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start costnloss2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Costnloss2
%this is a modified version of Costnloss to determine the cost and losses
%for buildable designs - thus "2" subscripts
%the method remains much the same, accept the n2 matrix is passed in to the function
%and it returns just the cost2 and actual_loss2 vectors.
% Returns vectors of costs, losses, and stranding for each of the calculated
% designs where each row designates a different design.
% Inputs are
% 'Irms' (root mean squared current)
% 'lt' (turn length)
% 'N' (number of turns)
% 'awg' (a vector of strand diameters for each of the optimized designs)
% 'rhoc' (the resistivity of copper at operating temperature)
% 'kl' (a parameter defined by time-rate of change of currents and the modified dynamic loss matrix).
% The outputs
% cost2 - are an array called 'cost' (Cost function of litz wire is assumed
% to be proportional to '(1+(k/(DC^4))+(kk))*n*lt*N' where 'k' is 1.1e-26,
% 'kk' is 2e-9, 'n' is the number of strands in the design in question
% 'lt' is the length of a single turn, and 'N' is the number of turns.
% where, dc is strand diameter, n is strand number, lt is turn
% length, and N is turn number)
%and
% 'actual_loss2' (the total winding power lost in each design and is equal to 'Irms^2*lt*N*rhoc*Fr/n/(pi*dc^2/4)' where all
% parameters are as defined above and rhoc is the resistivity of copper at
% operating temperature), and 'n' is defined to be the number of strands implemented
% in each design.
function [cost2,actual_loss2]= costnloss2(Irms,lt,n2,N,awg,rhoc,kl)
%we need to calculate the optimal values for Fr at different guages, so do that
%and put result into Fr[]
%Fr=1.+1./(1.-((2*COSTm(dc))./(dc.*dCOSTm(dc)))); % gives optimal Fr but we don't want to calculate optimal
% desgins anymore.
%n = numstrands(dc,kl,Fr); %numstrands will solve equation (2)
%now n holds a vector of the different number of strands for optimal Fr at
%different wire guages
%now calculate actual_loss
[numwind,numdesigns] = size(awg);
dc = awgtomet(awg); %convert wire guage to meters for calculation
Fr = [];
actual_loss = [];
klm = kl*ones(1, numdesigns);
%find Fr
Fr = 1 +klm.*n2.^2.*dc.^6.*(pi/4).^3;
%calculate the loss
actual_loss2 = (((Irms.^2).*lt.*N.*rhoc)*ones(1, numdesigns)).*Fr./n2./(pi.*dc.^2/4);
if length(lt) == 2
%this will print out the dc resistnace of the first two windings.
% Rdc1 = (lt(1).*N(1).*rhoc)./(n2(1).*(pi.*dc(1).^2/4))
% Rdc2 = (lt(2).*N(2).*rhoc)./(n2(2).*(pi.*dc(2).^2/4))
end
% now calculate cost for the cost/loss plot
Co = 0;
cm = COSTm(dc);
cost2 = (Co+(cm.*(dc.^2).*n2)).*( (lt.*N)*ones(1,numdesigns) );
%%%%%%%End of Costnloss2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start of fullbobbindesigns%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [n2] = fullbobbindesigns(awg, a, n, fp, N, userfp, bb, insulbuild)
%this function will determine the full-bobbin designs using the insulation build
%the inputs:
% awg - a vector of awg wires sizes
% n - a vector of the number of strands in each design
% fp - the packing factor required to build an optimal design
% N - a vector of the number of turns in each winding
% userfp - the packing factor specified by the user
% bb - the breadth of the bobbin window
% a - a vector containing the height of each winding
% insulbuild - the insulation build selected by user, either single or
% heavy
numwind = length(a);
numdesigns=length(awg);
if insulbuild == 's'
diameter = singlebuild(awg);
else
diameter = heavybuild(awg);
end
for i=1:numwind
for j=1:numdesigns
%if fp(i,j)>userfp
n2(i,j) = userfp*bb*a(i)/((diameter(j)^2)*N(i));
%else
% n2(i,j) = n(i,j);
%end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end fullbobbindesigns%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55start of singlebuild%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dt] = singlebuild(awg)
% This function determines the total diameter (in meters) of wire with single build insulation
% using the relationship determined by Charles Sullivan in "Optimal Choice for Number
% of Strands in a Litz Wire TRansformer Winding" published in IEEE Transaction on Power
% Electronics, Vol 14, No. 2, March 1999
%for wire sizes between 30-60 AWG
%the inputs:
% awg - a vector of awg wire sizes
%the output:
% dt - a vector of the total diameter of wire with insulation
alpha = 1.12; %constant for equation of single build insulation with AWG 40 as reference
beta = .97; %constant for equation of single build insulation with AWG 40 as reference
dr = .000079; %the reference diameter for function of AWG 40 in meters
dc = awgtomet(awg);
dt = dr*alpha*((dc./dr).^beta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end singlebuild%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start heavybuild%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dt = heavybuild(awg)
% This function determines the total diameter of wire with single build insulation
% using the relationship determined by Charles Sullivan in "Optimal Choice for Number
% of Strands in a Litz Wire TRansformer Winding" published in IEEE Transaction on Power
% Electronics, Vol 14, No. 2, March 1999
%the inputs:
% awg - a vector of awg wire sizes
%the output:
% dt - a vector of the total diameter of wire with insulation
alpha = 1.24; %constant for equation of single build insulation with AWG 40 as reference
beta = .94; %constant for equation of single build insulation with AWG 40 as reference
dr = .000079; %the reference diameter for function of AWG 40 in meters
dc = awgtomet(awg);
dt = dr*alpha*((dc./dr).^beta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end heavybuild%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5555
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start smartround%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r] = smartround(number,sigdigits)
%this function allows a number or vector to be rounded to a specified number of significant digits
%the inputs:
% number - a vector of the numbers you wnat to round
% sigdigits - the number of significant digits you want
%the output:
% r- a vector of rounded numbers
r = [];
[row, col] = size(number);
for i = 1:row
for j = 1:col
leadingdigit =floor(log10(number(i,j)))+1;
n = leadingdigit - sigdigits;
r(i,j) = 10.^n.*round(number(i,j).*10.^-n);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end smartround%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start optdesign table%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function htmltable = optdestable(optdesignlabels, totcost, totloss, numwind, n, dia, fp)
%generate table for optimal design and their packing factors
%the inputs:
% dia - a vector of of wire diameters
% n - a vector of the number of strands in each design
% fp - the packing factor required to build an optimal design
% totcost - the total cost of the optimal designs
% totloss - the total loss of the optimal designs
% optdesignlabels - the labels for the optimal design
title = {'Optimal Design Table'};
disp(title)
totcost = smartround(totcost, 3);
totloss = smartround(totloss, 3);
totcost = totcost';
totloss = totloss';
gaugeheading = {'Gauge','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge'; ' (W1)', ' (W2)',' (W3)',' (W4)',' (W5)',' (W6)', ' (W7)',' (W8)',' (W9)',' (W10)'};
gaugeheading = gaugeheading(:, 1:numwind);
heading = {'Relative','Loss','NumStrands',' NumStrands','NumStrands','NumStrands','NumStrands','NumStrands',' NumStrands','NumStrands','NumStrands','NumStrands'; 'Cost', 'in Watts', ' (W1)', ' (W2)',' (W3)',' (W4)',' (W5)',' (W6)', ' (W7)',' (W8)',' (W9)',' (W10)'};
heading = heading(:, 1:(2+numwind));
heading2 = {'Packing','Packing','Packing','Packing','Packing','Packing','Packing','Packing','Packing','Packing'; 'Factor (W1)','Factor (W2)','Factor (W3)','Factor (W4)','Factor (W5)','Factor (W6)','Factor (W7)','Factor (W8)','Factor (W9)','Factor (W10)'};
heading2 = heading2(:, 1:numwind);
heading = cat(2,gaugeheading, heading);
%heading = cat(2,heading, heading2);
n = smartround(n,3);
fp = smartround(fp,3);
n = n';
dia = dia';
fp = fp';
htmltable = [dia totcost totloss n];
htmltable = num2cell(htmltable);
%optdesignlabels = cellstr(optdesignlabels);
%htmltable = [optdesignlabels, htmltable];
htmltable = [heading; htmltable];
s.htmltable = htmltable; %#ok<*STRNU>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end optdesigntable%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start buildabletable%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bdtable = buildabletable(bbdesignlabels, totcost2round, totloss2round, numwind,n2round,awground)
%generate Html table for buildable design %uses the smartround function to round cost and loss columns
%the inputs:
% awground - the rounded awg wire sizes of buildable deigns
% n2round - a vector of the number of strands in each buildable design
% totcost2round - the total cost of the buildable designs
% totloss2round - the total loss of the buildable designs
% bbdesignlabels - the labels for the buildable designs
%create heading for the table
bdtitle = {'Buildable Design Table'};
disp(bdtitle)
bdcost = totcost2round;
bdcost = smartround(bdcost,3);
bdcost = bdcost';
gaugeheading = {'Design','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge','Gauge'; 'Number','(W1)', '(W2)','(W3)','(W4)','(W5)','(W6)', '(W7)','(W8)','(W9)','(W10)'};
gaugeheading = gaugeheading(:, 1:1+numwind);
bdheading = {'Relative','Loss','NumStrands',' NumStrands','NumStrands','NumStrands','NumStrands','NumStrands',' NumStrands','NumStrands','NumStrands','NumStrands'; 'Cost', 'in Watts', '(W1)', '(W2)','(W3)','(W4)','(W5)','(W6)', '(W7)','(W8)','(W9)','(W10)'};
bdheading = bdheading(:, 1:(2+numwind));
bdheading = cat(2,gaugeheading,bdheading);
%get data in porper format for the table
n2round = n2round';
[totloss2round] = smartround(totloss2round, 5);
totloss2round = totloss2round';
awground = awground';
bdtable = [awground bdcost totloss2round n2round];
bdtable = num2cell(bdtable);
bbdesignlabels = cellstr(bbdesignlabels);
bdtable = [bbdesignlabels, bdtable];
bdtable = [bdheading; bdtable];
s.bdtable = bdtable;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end buildabletable%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start roundesigns$%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [nrounded, awground] = rounddesigns(awg,n2, numwind) %#ok<*DEFNU>
%this program determines how to round designs and number of designs to create
%the inputs:
% numwind - the number oif windings in the component
% awg - a vector of awg sizes
% n2 - the number of strands in each buildable design
%the outputs:
% nrounded - the rounded vector of strandings
% awground - the vector of rounded awg sizes
numdesign = length(awg);
count = zeros(numwind,numdesign);
high = zeros(numwind,numdesign);
low = zeros(numwind,numdesign);
% determines the remainder of buildable designs then
% determines uses the remainder then rounds down if remaindern2 < 0.25
% rounds up if remaindern2>0.75 and creates both if remaindern2 is between
% 0.25 and 0.75
for k =1:numwind
for i=1:numdesign
if n2(k,i) <= 1
n2(k,i) = ceil(n2(k,i));
elseif n2(k,i) >= 10
n2(k,i) = round(n2(k,i));
else
remainder = n2(k,i) - floor(n2(k,i));
if remainder < 0.25
n2(k,i) = floor(n2(k,i));
elseif remainder <0.75
count(k,i) = 1;
low(k,i) = floor(n2(k,i));
high(k,i) = ceil(n2(k,i));
else
n2(k,i) =ceil(n2(k,i));
end
end
end
end
% this section of code creates the
% rounded n matrix and awg vectors
% totaldesigns = sum(count)+numdesign;
for k = 1:numwind
j = 1;
for i = 1:numdesign
if low(k,i) == 0 && high(k,i) == 0
awground(k,j) = awg(i);
nrounded(k,j) = n2(k,i);
j = j+1;
elseif low(k,i) ~= 0 || high(k,i) ~=0
awground(k,j) = awg(i);
awground(k,j+1) = awg(i);
nrounded(k,j) = low(k,i);
nrounded(k,j+1) = high(k,i);
j = j+2;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end%%rounddesigns%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%begin testnan%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [inputvector, kcostmatrix, klossmatrix] = testnan(inputvector, kcostmatrix, klossmatrix)
%This function prunes matrix that contains NaN that result from trying to extrapolate
%from the lookup table
%the inputs:
% inputvector -
% kcostmatirx - the cost matrix that you want to remove Nan from
% klossmatrix - the loss matrix that you want to remove Nan from
%the outputs:
% inputvector -
% kcostmatirx - the cost matrix without any NaN
% klossmatrix - the loss matrix without any NaN
nanvector = isnan(inputvector);
[row, columns] = size(nanvector);
for z = 1:columns
for z2 = 1:row
if nanvector(z2,z) == 1
inputvector(:,z) = [];
kcostmatrix(:,z) = [];
klossmatrix(:,z) = [];
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end testnan %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%begin testnanloss%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [inputvector, costW, kcostmatrix, klossmatrix] = testnanloss(inputvector, costW, kcostmatrix, klossmatrix)
%This function prunes matrix that contains NaN that result from trying to extrapolate
%from the lookup table
%the inputs:
% inputvector -
% kcostmatirx - the cost matrix that you want to remove Nan from
% klossmatrix - the loss matrix that you want to remove Nan from
% costW - the cost matrix that you want to remove Nan from
%the outputs:
% inputvector -
% kcostmatirx - the cost matrix without any NaN
% klossmatrix - the loss matrix without any NaN
% costW - the cost matrix without any Nan
count = 0;
nanvector = isnan(inputvector);
[row, columns] = size(nanvector);
for z = 1:columns
for z2 = 1:row
if nanvector(z2,z) == 1
inputvector(:,z-count) = [];
costW(:,z-count) = [];
kcostmatrix(:,z-count) = [];
klossmatrix(:,z-count) = [];
count = count + 1;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end testnanloss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%begin testnandia%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [inputvector, costW, lossW, kcostmatrix, klossmatrix] = testnandia(inputvector, costW, lossW, kcostmatrix, klossmatrix)
%This function prunes matrix that contains NaN that result from trying to extrapolate
%from the lookup table
%the inputs:
% inputvector -
% kcostmatirx - the cost matrix that you want to remove Nan from
% klossmatrix - the loss matrix that you want to remove Nan from
% costW - the cost matrix that you want to remove Nan from
% lossW - the loss matrix that you want to remove Nan from
%the outputs:
% inputvector -
% kcostmatirx - the cost matrix without any NaN
% klossmatrix - the loss matrix without any NaN
% costW - the cost matrix matrix without any NaN
% lossW - the loss matrix matrix without any NaN
count = 0;
nanvector = isnan(inputvector);
[row, columns] = size(nanvector);
for z = 1:columns
for z2 = 1:row
if nanvector(z2,z) == 1
inputvector(:,z-count) = [];
costW(:,z-count) = [];
lossW(:,z-count) = [];
kcostmatrix(:,z-count) = [];
klossmatrix(:,z-count) = [];
count = count + 1;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end testnandia %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start roundtoAWGrange%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function roundedAWG = roundtoAWGrange(dia, awg)
% this function rounded the matrix of wire diameters to the nearest value in awg
% vector passed in
% created by JDP 7/29/03
%the inputs:
% dia - a vector of diamters in meters
% awg - the vector of awg sizes you want to round to
%the output:
% roundedAWG - a vector of rounded AWG sizes that are closest to the
% diameters passed in
%determine size of dia matrix and setup roundedAWG matrix
[numwind, numdesigns] = size(dia);
roundedAWG = zeros(numwind, numdesigns);
%loop thru dia matrix and round nearest value in awg vector
for i = 1:numwind
for j=1:numdesigns
[minvalue, index] = min(abs(dia(i,j)-awg));
roundedAWG(i,j) = awg(index);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end of roundtoAWGrange%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start originaldesign%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tot_orig_cost, tot_orig_loss] = originaldesign(orig_numstrands, orig_wiresize, Irms, len, N, rhoc, kl, cost44)
% find cost and loss for original component design
% returns total origal cost and total original loss
%the inputs:
% Irms - a vector of the rms value of the current in each winding
% len - the average turn length for each winding
% bw - the breadth of the winding window
% N - a vetor of the number of turns in each winding
% numwind - the number of windings in the component
% rhoc - the resistivity of copper at given temperature
% kl - a parameter defined by time-rate of change of currents and the modified dynamic loss matrix
% cost44 - the cost used to mormalize the designs to awg 44
if orig_numstrands ~= 0
% Irms
% len
% orig_numstrands
% N
% orig_wiresize
% kl
[orig_cost, orig_loss]= costnloss2(Irms, len, orig_numstrands, N, orig_wiresize, rhoc, kl);
tot_orig_cost = sum(orig_cost)/cost44;
tot_orig_loss = sum(orig_loss);
else
tot_orig_cost = 0;
tot_orig_loss = 0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end orignaldesign%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%start windinght1d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [a, hb, bb, h] = windinght1d(Irms, bw, areabb, N, numwind)
%this function is designed to determine a, hb, bb and h for the 1d versions
%the inputs:
% Irms - a vector of the rms value of the current in each winding
% bw - the breadth of the winding window
% N - a vetor of the number of turns in each winding
% numwind - the number of windings in the component
%the outputs:
% a - Width of Winding Cross Section for each winding (in h direction)
% hb - the height of the bobbin
% bb - the breadth of the bobbin
% h - the height of the core window
%this is done by : the height of the winding = ((number of turns)(Irms)/(sum of N*Irms for allwindings))*height of bobbin
%function starts here
bb = bw;
hb = areabb/bw;
h = hb;
sumNIrms= sum(N.*Irms);
%determine ht of each winding
for i=1:numwind
hwinding(i) = (N(i)*Irms(i)/sumNIrms)*hb;
end
a = hwinding;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end windinght1d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%