%
% PriusEfficiency.m
%
% Prius HSD efficiency charts in Octave / Matlab. It all started with the question "how to minimize the flow over the electric path"
% (See http://www.priusfreunde.de/portal/index.php?option=com_kunena&Itemid=117&func=view&catid=13&id=236212 [German])
%
% The basic concept is that, assuming battery current is zero, all power going over the electic path is either generated by MG1
% consumed by MG2 or vice versa. The mechanical power of MG1 can be calculated by multiplying RPM and torque. MG1 RPM can be calculated
% from ICE RPM and vehicle speed (MG2 RPM is linear speed). MG1 torque can be calculated from ICE torque: MG1 torque is always ~28% of
% the ICE torque. Assuming that the ICE follows its "basic operating line" we know the ICE torque for each ICE RPM. This is all we
% need to calculate for combination of ICE RPM and vehicle speed the ICE power and MG1 power. Comparing those two values gives us an
% idea how the power is distributed between the mechanical and electric link.
%
% Licence CC BY-NC-SA 4.0, author proprius, contact me at http://www.priusfreunde.de or http://priuschat.com
%
% Version 1.0, 2014-02-02, proprius, Initial version for NHW20
% 
% The script makes some assumptions and simplifications. So reality will be a bit different. Here is what I can think of:
%
% - Battery current is assumed to be zero
% - Electric path efficiency (current between MG1 and MG2 over the inverter) is assumed to be 80%. This could be made more
%   precice by looking up the efficiency based on e. g. torque, RPM, current, and (boost) voltage. (see OakRidge_2004Prius.pdf for a start)
% - Maximum torque for MG1 when running as generator is unknown, especially at very low RPM. In theory (with 100% efficiency) it would not
%   require any power to generate a torque at 0 RPM. In reality there will be some current required and all of it will will be converted to
%   heat in MG1. The question is: At which speed is the net current to and from MG1 zero and what is the torque in that mode?
%   My very limited test driving suggests that this might be around 800 RPM, because that is when MG2 torque is at an absolute minimum, suggesting
%   that there is no flow over the electric path. Please test by watching MG1 RPM, MG2 torque, ICE RPM and battery current with your favorite
%   CAN bus tool (I have X-Gauges for ScanGauge).
% - Other values like rolling and air resistance, fuel and air density, and weight are set to default values. In reality they depend
%   on many factors, e. g. temperature, tires, etc.
% - Constant speed plot is for flat ground and no wind
% - Losses from gears and chain not included
% - Torque/BSFC values are taken from chart in Toyota paper SAE 2004-01-0064,
%   "Development of New-Generation Hybrid System THS II - Drastic Improvement of Power Performance and Fuel Economy,"
%   (see http://priuschat.com/threads/bsfc-challenge.55947/ )
%   + It is assumed that the basic operating line of followed, which in reality is not always true, see Fig. 22 in the link above.
%   + ICE is assumed to be warmed up
%   + ICE is assumed to be >= 1200 RPM and ICE power > 9.4 kW all the time
%   + The original chart only goes up to 4500 RPM. I added values for 5000 RPM based on the max output power of 57 kW
%     for the ICE. The BSCF at 5000 RPM is just a good guess based on the 4500 RPM chart.
%
% Resources (Thanks to all the authors):
%
% http://eahart.com/prius/psd/
% http://prius.ecrostech.com/original/Understanding/PowerSplitDevice.htm
% http://priuschat.com/threads/bsfc-challenge.55947/
% http://home.arcor.de/sg_temp/prius/fahrprofilrechner.htm
% http://priuschat.com/threads/mg1-rpm.59608/
% http://priuschat.com/threads/ice-to-wheel-transmission-efficiency-value-at-60-mph.115094/
% http://priuschat.com/threads/introduction-to-prius-power-flow.29352/
% http://priuschat.com/threads/two-mg1-questions.117253/%post-1712003
% http://ecee.colorado.edu/~ecen5017/notes/OakRidge_2004Prius.pdf
% http://info.ornl.gov/sites/publications/files/Pub26762.pdf

global MODEL ELECTRIC_PATH_EFFICIENCY SUN_GEAR_TEETH RING_GEAR_TEETH SUN_TO_RING SUN_TO_CARRIER MG1_ICE_TORQUE_FACTOR;
global MAX_RPM_MG1 MIN_RPM_MG1 MAX_RPM_MG2 MAX_KMH_FOR_MG2 MAX_KMH_SPEEDOMETER SPEEDOMETER_FACTOR MIN_RPM_ICE MAX_RPM_ICE;
global ICE_RPM_TABLE ICE_TORQUE_TABLE ICE_BSFC_TABLE FUEL_DENSITY;
global CW FRONTAREA ROLLING WEIGHT AIR_DENSITY;

MODEL                    = 'Toyota Prius 2004-2009 (NHW20)';
ELECTRIC_PATH_EFFICIENCY = 0.80;
SUN_GEAR_TEETH           = 30;
RING_GEAR_TEETH          = 78;
SUN_TO_RING              = RING_GEAR_TEETH ./ SUN_GEAR_TEETH;                           % 2.6
SUN_TO_CARRIER           = SUN_TO_RING .+ 1;                                            % 3.6
MG1_ICE_TORQUE_FACTOR    = -(SUN_GEAR_TEETH ./ (SUN_GEAR_TEETH .+ RING_GEAR_TEETH));    % ~ -0.2777
MAX_RPM_MG1              = 10000;
MIN_RPM_MG1              = -MAX_RPM_MG1;
MAX_RPM_MG2              = 6400;
MAX_KMH_FOR_MG2          = 181.9;
SPEEDOMETER_FACTOR       = 1.075;
MAX_KMH_SPEEDOMETER      = 170 .* SPEEDOMETER_FACTOR;
MIN_RPM_ICE              = 1200;
MAX_RPM_ICE              = 5000;
ICE_RPM_TABLE            = [1200 1290 1700 2000 2220   2800 3230 3500 4500 5000];
ICE_TORQUE_TABLE         = [  75   77   81   85   86.5   91   94   97   97  109];
ICE_BSFC_TABLE           = [ 242  240  235  232  230    230  230  230  235  240];
FUEL_DENSITY             = 0.750;
CW                       = 0.26;
FRONTAREA                = 2.5;
ROLLING                  = 0.01;
WEIGHT                   = 1375;
AIR_DENSITY              = 1.225;

function retval = gramToLiter (gram)
  global FUEL_DENSITY;
  retval = gram ./ (FUEL_DENSITY .* 1000);
end  

function retval = kmhToMeterPerSec (kmh)
  retval = kmh ./ 3.6;
end

function retval = speedometerToTrueSpeed (kmh)
  global SPEEDOMETER_FACTOR;
  retval = kmh ./ SPEEDOMETER_FACTOR;
end

function retval = powerConstantSpeedOnFlatGround (kmh)
  % Derived from (c) Hjermind 2008-2009 http://home.arcor.de/sg_temp/prius/fahrprofilrechner.htm
  global CW FRONTAREA ROLLING WEIGHT AIR_DENSITY;
  v = kmhToMeterPerSec(kmh);
  air = v .* (AIR_DENSITY ./ 2) .* CW .* FRONTAREA .* v .* v ./ 1000;
  roll = v .* ROLLING .* WEIGHT .* 9.81 ./ 1000;
  retval = air .+ roll; 
end

function retval = kmhToRpmMg2 (kmh)
  global MAX_RPM_MG2 MAX_KMH_FOR_MG2; 
  retval = (kmh .* MAX_RPM_MG2) ./ MAX_KMH_FOR_MG2;
end  

function retval = rpmMg2ToKmh (rpmMg2)
  global MAX_RPM_MG2 MAX_KMH_FOR_MG2;
  retval = MAX_KMH_FOR_MG2 .* (rpmMg2 ./ MAX_RPM_MG2);
end  

function retval = rpmMg1 (rpmMg2,rpmIce)
  global SUN_TO_RING SUN_TO_CARRIER;
  retval = (SUN_TO_CARRIER .* rpmIce) .- (SUN_TO_RING .* rpmMg2);
end

function retval = rpmMg2 (rpmMg1,rpmIce)
  global SUN_TO_RING SUN_TO_CARRIER;
  retval = ((SUN_TO_CARRIER .* rpmIce) .- rpmMg1) ./ SUN_TO_RING;
end

function retval = rpmMg1IceToKmh (rpmMg1,rpmIce)
  retval = rpmMg2ToKmh(rpmMg2(rpmMg1,rpmIce));
end

function retval = powerMg1relativeToIce (rpmMg2,rpmIce)
  global MG1_ICE_TORQUE_FACTOR;
  retval = (rpmMg1(rpmMg2,rpmIce) ./ rpmIce) .* MG1_ICE_TORQUE_FACTOR;
end

function retval = powerMg1relativeToIceKmh (kmh, rpmIce)
  retval = powerMg1relativeToIce(kmhToRpmMg2(kmh),rpmIce);
end  

function [torque bsfc] = torqueAndBsfcIce (rpmIce)
  global ICE_RPM_TABLE ICE_TORQUE_TABLE ICE_BSFC_TABLE;
  torque = interp1(ICE_RPM_TABLE,ICE_TORQUE_TABLE,rpmIce);
  bsfc = interp1(ICE_RPM_TABLE,ICE_BSFC_TABLE,rpmIce);
end

function torque = torqueIce (rpmIce)
  [torque bsfc] = torqueAndBsfcIce(rpmIce);
end

function bsfc = bsfcIce (rpmIce)
  [torque bsfc] = torqueAndBsfcIce(rpmIce);
end

function retval = powerIce (rpmIce)
  retval = (torqueIce(rpmIce) .* rpmIce) ./ 9549.3;
end

function retval = powerMg1 (kmh,rpmIce)
  retval = powerIce(rpmIce) .* powerMg1relativeToIceKmh(kmh,rpmIce);
  
  % When regenerated power is too low (because MG1 rpm is too low), then power is consumed to slow down MG1
  minPower = -2;
  regenerationPowerTooLow = (retval >= minPower && retval <= 0);
  retval = retval .* (1 .- regenerationPowerTooLow) .+ minPower .* regenerationPowerTooLow;
end  

function retval = powerAfterPsd (kmh,rpmIce)
  global ELECTRIC_PATH_EFFICIENCY;
  retval = powerIce(rpmIce) .- (abs(powerMg1(kmh,rpmIce)) * (1 .- ELECTRIC_PATH_EFFICIENCY));
end 

function retval = powerSurplus (kmh,rpmIce)
  retval = powerAfterPsd(kmh,rpmIce) .- powerConstantSpeedOnFlatGround(kmh);
end 

function literPerHundredKm = consumption (kmh,rpmIce)
  literPerHour = gramToLiter(bsfcIce(rpmIce) .* powerIce(rpmIce));
  literPerHundredKm = (literPerHour .* 100) ./ kmh;
end

function retval = rpmMg1kmh (kmh, rpmIce)
  retval = rpmMg1(kmhToRpmMg2(kmh),rpmIce);
end  

function retval = minKmhSpeedometer (rpmIce)
  global MAX_RPM_MG1;
  retval = rpmMg1IceToKmh(MAX_RPM_MG1,rpmIce);
  retval = retval .* (retval >= 0);
end  

function retval = maxKmhSpeedometer (rpmIce)
  global MIN_RPM_MG1 MAX_KMH_SPEEDOMETER;
  retval = rpmMg1IceToKmh(MIN_RPM_MG1,rpmIce);
  retval = retval .* (retval <= MAX_KMH_SPEEDOMETER) + MAX_KMH_SPEEDOMETER .* (retval > MAX_KMH_SPEEDOMETER);
end  

function retval = efficiency (kmh,rpmIce)
  global ICE_BSFC_TABLE;
  bsfc = bsfcIce(rpmIce);
  retval = 100 .* (powerAfterPsd(kmh,rpmIce) ./ powerIce(rpmIce)) ./ (bsfc ./ min(ICE_BSFC_TABLE));
end  

function plotConstantSpeed (X,Y,lineColor,textColor)
  Z = powerSurplus(speedometerToTrueSpeed(X),Y);
  [C,h] = contour(X,Y,Z,[0,0]);
  clabel(C,h,'rotation',45);
  line = findobj(h,'String','0');
  set(line,'Color',textColor); 
  set(line,'String','Constant Speed');
  set(h,'LineWidth',3);
  set(h,'LineColor',lineColor);
end

function plotRpmMg1 (X,Y,rpm,lineColor,textColor)
  Z = rpmMg1kmh(speedometerToTrueSpeed(X),Y);
  [C,h] = contour(X,Y,Z,[rpm,rpm]);
  clabel(C,h,'rotation',45,'LabelSpacing',170);
  line = findobj(h,'String',sprintf('%d',rpm));
  set(line,'Color',textColor);
  set(line,'String',sprintf('MG1 %d rpm',rpm));
  set(h,'LineWidth',3);
  set(h,'LineColor',lineColor);
end  

function plotMinus2kwMg1 (X,Y)
  Z = powerMg1(speedometerToTrueSpeed(X),Y);                 
  [C,h] = contour(X,Y,Z,[-2,-2]);
  clabel(C,h,'rotation',45);
  line = findobj(h,'String','-2');
  set(line,'Color','w'); 
  set(line,'String','MG1 -2 kW');
  set(h,'LineWidth',3);
end

function axisSpeedometerAndRpm ()
  global MAX_KMH_SPEEDOMETER MIN_RPM_ICE MAX_RPM_ICE;
  xlabel('Speed as displayed on dashboard (km/h)');
  set(gca,'XTick',0:10:MAX_KMH_SPEEDOMETER);
  ylabel('ICE RPM');
  set(gca,'YTick',[MIN_RPM_ICE MAX_RPM_ICE:-500:MIN_RPM_ICE]);
end

function licence (x,y)
  hLicence = text(x,y,'Created by proprius (CC BY-NC-SA 4.0)','clipping','off');
  set(hLicence,'FontSize',6);
end

function plotEfficiency (X,Y)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Efficiency map for ICE and electric path MG1<->Inverter<->MG2, 100 = no electric losses and', 'ICE most efficient. See accompanying text for assumptions and simplifications made.'});
  Z = efficiency(speedometerToTrueSpeed(X),Y);
  [C,h] = contourf(X,Y,Z,[50:100]);
  ch = clabel(C,h,'rotation',0);
  set(ch,'BackgroundColor',[1 1 .6])
  plotConstantSpeed(X,Y,'b','w');
  plotRpmMg1(X,Y,0,'g','w');
  plotRpmMg1(X,Y,800,[1.0 0.9 0.0],'w');
  licence(0,1020);
  hold off;
  set(gcf,'PaperUnits','inches','PaperPosition',[0 0 9 7]);
  saveas(hFigure,'PriusEfficiency.png');
end

function plotConsumption (X,Y)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Fuel consumption (L/100 km). See accompanying text for assumptions and simplifications made.'});
  Z = consumption(speedometerToTrueSpeed(X),Y);
  [C,h] = contourf(X,Y,Z,[0:30]);
  clabel(C,h,'rotation',0,'color','w');
  plotConstantSpeed(X,Y,'r','w');
  plotRpmMg1(X,Y,0,'g','w');
  plotRpmMg1(X,Y,800,[1.0 0.9 0.0],'w');
  licence(0,1020);
  hold off;
  set(gcf,'PaperUnits','inches','PaperPosition',[0 0 9 7]);
  saveas(hFigure,'PriusConsumption.png');
end

function plotPowerAfterPsd (X,Y)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Power after PSD (kW). See accompanying text for assumptions and simplifications made.'});
  Z = powerAfterPsd(speedometerToTrueSpeed(X),Y);
  [C,h] = contourf(X,Y,Z,[0:2:57]);
  clabel(C,h,'rotation',0);
  plotConstantSpeed(X,Y,'r','w');
  plotRpmMg1(X,Y,0,'g','w');
  plotRpmMg1(X,Y,800,[1.0 0.9 0.0],'w');
  licence(0,1020);
  hold off;
  set(gcf,'PaperUnits','inches','PaperPosition',[0 0 9 7]);
  saveas(hFigure,'PriusPowerAfterPsd.png');
end

function plotPowerSurplus (X,Y)
  global MODEL;
  hFigure = figure('Position',[0,0,1600,1200]);
  hold on;
  grid on;
  axisSpeedometerAndRpm();
  title({MODEL, 'Power surplus: Power after PSD minus power required for constant speed on flat ground (kW).', 'See accompanying text for assumptions and simplifications made.'});
  Z = powerSurplus(speedometerToTrueSpeed(X),Y);
  [C,h] = contourf(X,Y,Z,[-40:2:40]);
  clabel(C,h,'rotation',0);
  plotConstantSpeed(X,Y,'b','k');
  plotRpmMg1(X,Y,0,[0.0 0.5 0.0],'k');
  plotRpmMg1(X,Y,800,'k','k');
  licence(0,1020);
  hold off;
  set(gcf,'PaperUnits','inches','PaperPosition',[0 0 9 7]);
  saveas(hFigure,'PriusPowerSurplus.png');
end

function plotTorqueAndPower (rpmIceRange)
  global MODEL MIN_RPM_ICE MAX_RPM_ICE ICE_TORQUE_TABLE;
  hFigure = figure;
  hold on;
  grid on;
  title(MODEL);
  [hAxes,hLine1,hLine2] = plotyy(rpmIceRange,torqueIce(rpmIceRange),rpmIceRange,powerIce(rpmIceRange),'plot');
  set(hLine1,'LineWidth',2);
  set(hLine2,'LineWidth',2);
  set(hAxes(1),'XTick',[MIN_RPM_ICE MAX_RPM_ICE:-500:MIN_RPM_ICE]);
  set(hAxes(2),'XTick',0);
  set(hAxes(1),'YLim', [60 120]);
  set(hAxes(2),'YLim', [0 60]);
  set(hAxes(1),'YTick',[60:5:120]);
  set(hAxes(2),'YTick',[0:5:60]);
  set(get(hAxes(1),'Xlabel'),'String','ICE RPM');
  set(get(hAxes(1),'Ylabel'),'String','ICE Torque (Nm)');
  set(get(hAxes(2),'Ylabel'),'String','ICE Power (kW)');
  
  colorTorque = 'b';
  colorPower = 'r';
  set(hAxes(1),'ycolor',colorTorque); % not working
  set(get(hAxes(1),'YLabel'),'color', colorTorque); 
  set(hLine1,'Color',colorTorque);
  set(hAxes(2),'ycolor',colorPower); % not working
  set(get(hAxes(2),'YLabel'),'color', colorPower); 
  set(hLine2,'Color',colorPower);
  
  hold off;
  set(gcf,'PaperUnits','inches','PaperPosition',[0 0 5.5 4]);
  saveas(hFigure,'PriusTorqueAndPower.png');
end
  
rpmIceRange = MIN_RPM_ICE:.2:MAX_RPM_ICE; kmhSamples = 1000;  % very high quality plot
%rpmIceRange = MIN_RPM_ICE:20:MAX_RPM_ICE; kmhSamples = 100;   % mid quality plot
%rpmIceRange = MIN_RPM_ICE:400:MAX_RPM_ICE; kmhSamples = 5;    % low quality plot

X = linspace(minKmhSpeedometer(rpmIceRange),maxKmhSpeedometer(rpmIceRange),kmhSamples);
Y = repmat((rpmIceRange)',1,kmhSamples);

%Z = powerMg1(speedometerToTrueSpeed(X),Y);                  [C,h] = contourf(X,Y,Z,[-30:2:30]);                     clabel(C,h,'rotation',0);
%Z = powerMg1relativeToIceKmh(speedometerToTrueSpeed(X),Y);  [C,h] = contourf(X,Y,Z,[0:0.1:2.5]);                    clabel(C,h,'rotation',0);
%Z = rpmMg1kmh(speedometerToTrueSpeed(X),Y);                 [C,h] = contour (X,Y,Z,[-MIN_RPM_MG1:400:MAX_RPM_MG1]); clabel(C,h,'rotation',0);

plotTorqueAndPower(rpmIceRange);
plotEfficiency(X,Y);
plotConsumption(X,Y);
%plotPowerAfterPsd(X,Y);
plotPowerSurplus(X,Y);

pause