%% wiebeFit.m % % [FitCoeffs,CAparams,Wiebes] = wiebeFitBF(NumWiebes,Breakdown,CAs,HRR,CAparams,Guesses,Ranges,VecSize,a) % % Fits a Wiebe function to a heat release rate curve via Brute Force % % NumWiebes = # of sumltaneously brute-forced Wiebe functions, 1, 2, or 3 % CAs: Vector of crank angles % HRR: Vector of heat release rates % CAparams: 3-element vector, 1st element = start of CA window, 2nd % element = end of CA window, 3rd element = CA interval % Guesses: NumWiebes-by-4 matrix of starting guesses for the wiebe % function parameters, 1st col = theta_0, 2nd col = delta theta, % 3rd col = Q, 4th col = M % VecSize = Number of possibilities to try for each parameter % EfficiencyFactor = Efficiency factor "a" in the Wiebe functions function [FitCoeffs,CAparams,Wiebe1s] = wiebeFit(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor) options = optimoptions('fmincon','FiniteDifferenceStepSize', 0.1); [BestCoeffs, fval] = fmincon(@(CoeffArray) OBJ(CoeffArray,NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor), [Guesses(:,1), Guesses(:,2), Guesses(:,3), Guesses(:,4)], options); % BestCoeffs(1,:) = cell2mat(CoeffArray(1,:)); %converts a cell array into an ordinary array (array in matrix) FitCoeffInds = fliplr(BestCoeffs); %returns BestCoeffs with its columns flipped in the left-right direction FitCoeffs = zeros(size(BestCoeffs)); for WiebeNum = 1:NumWiebes for CoeffNum = 1:4 FitCoeffs(WiebeNum,CoeffNum) = Coeffs(WiebeNum,CoeffNum,FitCoeffInds(WiebeNum,CoeffNum)); end end [Err, Wiebels] = OBJ(BestCoeffs, NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor); end function [Errs, Wiebels] = OBJ(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor) NewCAs = CAparams(1):CAparams(3):CAparams(2); NewHRR = zeros(1,length(NewCAs)); for CAnum = 1:length(NewCAs) CA = NewCAs(CAnum); [~,CAind] = min(abs(CAs-CA)); %locate the index Err = CAs(CAind)-CA; if Err == 0 NewHRR(CAnum) = HRR(CAind); else if CAind < 2 Delta = HRR(CAind+1)-HRR(CAind); InterpFac = Err/(CAs(CAind+1)-CAs(CAind)); %InterpolationFactor to calculate Err up to 1 deg else Delta = HRR(CAind)-HRR(CAind-1); InterpFac = Err/(CAs(CAind)-CAs(CAind-1)); end NewHRR(CAnum) = HRR(CAind)-Delta*InterpFac; end end NumCAs = length(NewCAs); NumCoeffs = NumWiebes*4; Ints = -1:(2/(VecSize-1)):1; %VecSize is Resolution = 11 Coeffs = []; Coeffs = ones(NumWiebes,4,VecSize); for SweepNum = 1:VecSize Coeffs(:,:,SweepNum) = Guesses+Ints(SweepNum)*Ranges; end Wiebes = calcWiebeFunction(NumWiebes,Coeffs,CAparams,EfficiencyFactor); % NumPoss = size(CoeffArray,1); NumPoss = 1; Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array if NumWiebes > 1; for Poss2 = 1:NumPoss Coeff2 = CoeffArray(Poss2,:); Wiebe2 = Wiebes(2,:,Coeff2{:}); NewWiebes = Wiebe1s+ones(NumPoss,1)*Wiebe2; for Poss3 = 1:NumPoss Coeff3 = CoeffArray(Poss3,:); Wiebe3 = Wiebes(3,:,Coeff3{:}); NewWiebes = NewWiebes+ones(NumPoss,1)*Wiebe3; Errs = sum((NewWiebes-ones(NumPoss,1)*NewHRR),2)/sum(NewHRR); end end else Errs = sum((Wiebe1s-ones(VecSize^4,1)*NewHRR).^2,2)/sum(NewHRR.^2); %sum(X,2) -> sum of each row [~,NewInds] = min(Errs); end end % function [FitCoeffs,CAparams,Wiebe1s] = wiebeFit(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor) % options = optimoptions('fmincon','FiniteDifferenceStepSize', 0.1); % [BestCoeffs, fval] = fmincon(@(CoeffArray) OBJ(CoeffArray,NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor), [Guesses(:,1), Guesses(:,2), Guesses(:,3), Guesses(:,4)], options); % % BestCoeffs(1,:) = cell2mat(CoeffArray(1,:)); %converts a cell array into an ordinary array (array in matrix) % % FitCoeffInds = fliplr(BestCoeffs); %returns BestCoeffs with its columns flipped in the left-right direction % FitCoeffs = zeros(size(BestCoeffs)); % for WiebeNum = 1:NumWiebes % for CoeffNum = 1:4 % FitCoeffs(WiebeNum,CoeffNum) = Coeffs(WiebeNum,CoeffNum,FitCoeffInds(WiebeNum,CoeffNum)); % end % end % [Err, Wiebels] = OBJ(BestCoeffs, NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor); % end % % function [Errs, Wiebels] = OBJ(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor) % NewCAs = CAparams(1):CAparams(3):CAparams(2); % NewHRR = zeros(1,length(NewCAs)); % for CAnum = 1:length(NewCAs) % CA = NewCAs(CAnum); % [~,CAind] = min(abs(CAs-CA)); %locate the index % Err = CAs(CAind)-CA; % if Err == 0 % NewHRR(CAnum) = HRR(CAind); % else % if CAind < 2 % Delta = HRR(CAind+1)-HRR(CAind); % InterpFac = Err/(CAs(CAind+1)-CAs(CAind)); %InterpolationFactor to calculate Err up to 1 deg % else % Delta = HRR(CAind)-HRR(CAind-1); % InterpFac = Err/(CAs(CAind)-CAs(CAind-1)); % end % NewHRR(CAnum) = HRR(CAind)-Delta*InterpFac; % end % end % NumCAs = length(NewCAs); % % NumCoeffs = NumWiebes*4; % Ints = -1:(2/(VecSize-1)):1; %VecSize is Resolution = 11 % Coeffs = []; % Coeffs = ones(NumWiebes,4,VecSize); % for SweepNum = 1:VecSize % Coeffs(:,:,SweepNum) = Guesses+Ints(SweepNum)*Ranges; % end % Wiebes = calcWiebeFunction(NumWiebes,Coeffs,CAparams,EfficiencyFactor); % % % NumPoss = size(CoeffArray,1); % NumPoss = 1; % % Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension % Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array % % Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension % Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array % % if NumWiebes > 1; % for Poss2 = 1:NumPoss % Coeff2 = CoeffArray(Poss2,:); % Wiebe2 = Wiebes(2,:,Coeff2{:}); % NewWiebes = Wiebe1s+ones(NumPoss,1)*Wiebe2; % for Poss3 = 1:NumPoss % Coeff3 = CoeffArray(Poss3,:); % Wiebe3 = Wiebes(3,:,Coeff3{:}); % NewWiebes = NewWiebes+ones(NumPoss,1)*Wiebe3; % Errs = sum((NewWiebes-ones(NumPoss,1)*NewHRR),2)/sum(NewHRR); % end % end % else % Errs = sum((Wiebe1s-ones(VecSize^4,1)*NewHRR).^2,2)/sum(NewHRR.^2); %sum(X,2) -> sum of each row % [~,NewInds] = min(Errs); % end % end % function [FitCoeffs,CAparams,Wiebe1s] = TwiebeFitBruteForce(NumWiebes,CAs,HRR,CAparams,Guesses,Ranges,VecSize,EfficiencyFactor) % % NewCAs = CAparams(1):CAparams(3):CAparams(2); % NewHRR = zeros(1,length(NewCAs)); % for CAnum = 1:length(NewCAs) % CA = NewCAs(CAnum); % [~,CAind] = min(abs(CAs-CA)); %locate the index % Err = CAs(CAind)-CA; % if Err == 0 % NewHRR(CAnum) = HRR(CAind); % else % if CAind < 2 % Delta = HRR(CAind+1)-HRR(CAind); % InterpFac = Err/(CAs(CAind+1)-CAs(CAind)); %InterpolationFactor to calculate Err up to 1 deg % else % Delta = HRR(CAind)-HRR(CAind-1); % InterpFac = Err/(CAs(CAind)-CAs(CAind-1)); % end % NewHRR(CAnum) = HRR(CAind)-Delta*InterpFac; % end % end % NumCAs = length(NewCAs); % % NumCoeffs = NumWiebes*4; % Ints = -1:(2/(VecSize-1)):1; %VecSize is Resolution = 11 % Coeffs = []; % Coeffs = ones(NumWiebes,4,VecSize); % for SweepNum = 1:VecSize % Coeffs(:,:,SweepNum) = Guesses+Ints(SweepNum)*Ranges; % end % Wiebes = calcWiebeFunction(NumWiebes,Coeffs,CAparams,EfficiencyFactor); % % CoeffArray = num2cell(permn(1:VecSize,4)); %converts permn... into an cell array 4 (11^4 different combinations) % NumPoss = size(CoeffArray,1); % % Wiebe1s = shiftdim(squeeze(Wiebes(1,:,:,:,:,:)),1); %shiftdim shifts the dimensions of squeeze... by 1 (to right side) - squeeze removes the singleton diemension % Wiebe1s = reshape(Wiebe1s,VecSize^4,NumCAs); %reshapes Wiebe1s into a VecSize^4-by-NumCAs array % % BestError = inf; % BestCoeffs = zeros(NumWiebes,4); % % if NumWiebes > 1; % for Poss2 = 1:NumPoss % Coeff2 = CoeffArray(Poss2,:); % Wiebe2 = Wiebes(2,:,Coeff2{:}); % NewWiebes = Wiebe1s+ones(NumPoss,1)*Wiebe2; % for Poss3 = 1:NumPoss % Coeff3 = CoeffArray(Poss3,:); % Wiebe3 = Wiebes(3,:,Coeff3{:}); % NewWiebes = NewWiebes+ones(NumPoss,1)*Wiebe3; % Errs = sum((NewWiebes-ones(NumPoss,1)*NewHRR),2)/sum(NewHRR); % [NewError,NewInds] = min(Errs); % if NewError < BestError % BestError = NewError; % BestCoeffs(1,:) = cell2mat(CoeffArray(NewInds,:)); %converts a cell array into an ordinary array % BestCoeffs(2,:) = cell2mat(CoeffArray(Poss2,:)); % BestCoeffs(3,:) = cell2mat(CoeffArray(Poss3,:)); % end % end % end % else % Errs = sum((Wiebe1s-ones(VecSize^4,1)*NewHRR).^2,2)/sum(NewHRR.^2); %sum(X,2) -> sum of each row % [~,NewInds] = min(Errs); % BestCoeffs(1,:) = cell2mat(CoeffArray(NewInds,:)); %converts a cell array into an ordinary array (array in matrix) % end % % FitCoeffInds = fliplr(BestCoeffs); %returns BestCoeffs with its columns flipped in the left-right direction % FitCoeffs = zeros(size(BestCoeffs)); % for WiebeNum = 1:NumWiebes % for CoeffNum = 1:4 % FitCoeffs(WiebeNum,CoeffNum) = Coeffs(WiebeNum,CoeffNum,FitCoeffInds(WiebeNum,CoeffNum)); % end % end % % end