load Tuning_data_Dec18 % original with a_dw=11.5 %load Tuning_data_alpha12 %load Tuning_data_alpha13 %xdata=[voltage_inputs(:,1),voltage_inputs(:,2),voltage_inputs(:,3),voltage_inputs(:,6),voltage_inputs(:,7),Mas_liqwater_anode]; %%% Try using anode section 1. %xdata=[y(:,[2 3 4 6 7 ]),ANODE_1(:,2)*V_P*VAPORMOLARMASS+ANODE_1(:,5)*V_P*WATERDENSITY]; % %%% Try using anode GDL total mass of water. %xdata=[y(:,[2 3 4 6 7 ]),... % ANODE_1(:,2)*V_P*VAPORMOLARMASS+ANODE_1(:,5)*V_P*WATERDENSITY+ANODE_2(:,2)*V_P*VAPORMOLARMASS+ANODE_2(:,5)*V_P*WATERDENSITY+ANODE_3(:,2)*V_P*VAPORMOLARMASS+ANODE_3(:,5)*V_P*WATERDENSITY]; %%% Try using Cathode section 1. %xdata=[y(:,[2 3 4 6 7 ]),CATHODE_1(:,3)*V_P*VAPORMOLARMASS+CATHODE_1(:,7)*V_P*WATERDENSITY]; %xdata=[y(:,[2 3 4 6 7 ]),CATHODE_1(:,7)*V_P*WATERDENSITY]; xdata=y(:,[2 3 4 6 7 8 ]); %ydata=interp1(time,AveCellV,tout,'linear')/1000; %ydata=interp1(time,sum(data(:,2:2+23),2)/24,tout,'linear')/1000; %Throw out the bad cell... with the cross leak %use everything %ydata=interp1(time,sum(data(:,2:2+23),2)/24,tout,'linear'); %Throw out the bad cell... with the cross leak % Finding index for min, max and median voltages, index selects time range % Dec 17 Data % index=find( y(:,7) > 1 ); % index=1:length(t); % %% Dec 18 Data index=find( t>300 & t< 3000 | t>8600& t< 15280 );%&y(:,7) > 1 % %% Dec 20 Data % index=find( ((t>3400 & t<13000) | (t>15302& t<22500 )| (t>13800& t<13900 ) )) ; %throw out the first bit of data % index=find( ((t>5200 & t<13000) | (t>15302& t<22500 ) )) ; %throw out the first bit of data %index=find( ((t>1000 & t<3240) |(t>3400 & t<13000) | (t>15302& t<22500 ) )) ; %index=find( (t>3400 & t<13000) ) ; % %% Dec 21 Data % index=find( ((t>4000 & t<16400) | (t>18500& t<24080 )| (t>24600& t<29350 ) )) ; %throw out the first bit of data % index=find( t>4000 & y(:,7)> 0 ) ; %throw out the first bit of data % index=find( ((t>4000 & t<13000) | (t>13440 & t<13590) | (t>14702& t<22500 ) )& y(:,7)>.1 ) ; %throw out the first bit of data % % index=find( ((t>4700 & t<13580) | (t>14700 & t<14950) | (t>15024& t<22500 ) )& y(:,7)>.1 ) ; %throw out the first bit of data % index=find( ((t>13800 & t<13900) |(t>3400 & t<13375) | (t>14700 & t<14950) | (t>15020& t<22500 )| (t>13460 & t<13580 ) | (t>29380 & t<55000 ) ) & y(:,7)>5) ; %throw out the first bit of data % index=find( ((t>3400 & t<13375) | (t>15048& t<22500 ) | (t>30000 & t<38700 )| (t>40700 & t<46200) | (t>47100 & t<51600) | (t>53100 & t<55300) ) & y(:,7)>5) ; %throw out the first bit of data % index=find( ((t>3400 & t<13375) | (t>15048& t<22500 ) | (t>30000 & t<38700 )| (t>40700 & t<46200) | (t>47100 & t<51600) ) & y(:,7)>5) ; %throw out the first bit of data %index=find( ((t>5180 & t<13375) | (t>15048& t<22500 ) | (t>30000 & t<37500 )| (t>40700 & t<46200) | (t>47100 & t<51600) ) & y(:,7)>5) ; %throw out the first bit of data size(index) Average_voltage=sum(data(:,2:2+23),2)/24; %Average Median_voltage=median(data(:,2:2+23),2); Error=data( : , 2:(2+23))-Median_voltage*ones(1,24); sqe=zeros(1,24); for i=1:24 % sqe(i)=(Error(:,i).')*Error(:,i); sqe(i)=(Error(index,i).')*Error(index,i); end [C,IMedian] = min(sqe,[],2); IMedian Median_voltage=data(:,1+IMedian); Min_voltage_LB=min( data( : , 2:(2+23) ).').'; Error=data( : , 2:(2+23))-Min_voltage_LB*ones(1,24); sqe=zeros(1,24); for i=1:24 sqe(i)=(Error(index,i).')*Error(index,i); end [C,IMin] = min(sqe,[],2); IMin Min_voltage=data(:,1+IMin); Max_voltage_UP=max( data( : , 2:(2+23) ).').'; Error=data( : , 2:(2+23))-Max_voltage_UP*ones(1,24); sqe=zeros(1,24); for i=1:24 sqe(i)=(Error(index,i).')*Error(index,i); end [C,IMax] = min(sqe,[],2); IMax Max_voltage=data(:,1+IMax); % vector of ur is formed from mscaledata %[PCaIn,TCaIn,RHCaIn,PAnIn,TAnIn,RHAnIn,... %1,2,3,4,5,6 % PCaOut,TCaOut,RHCaOut,PAnOut,TAnOut,RHAnOut,...%7,8,9,10,11,12 % WaCaIn,IFC,TFC,PurgeState]; %13,14,15,16 % find index for the end of the calibration data set LL=max(find(time= 1) = 1; AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - 2*AN_GDL_surf_blockage); %%%IGNORE TWL %Calculate the apparent current density using the stack current. i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] %E=-(-241980 -11.0132*(T_st-T0)+0.003968 *(T_st-T0).^2/2+1.12E-6*(T_st-T0).^3/3-T_st.*(-44.4 -11.0132*log(T_st./T0)+0.003968*(T_st-T0) +1.12E-6*(T_st-T0).^2/2) )/2/FARADAYS+R_G_U*T_st/2/FARADAYS.*log(p_H2.* sqrt(p_O2)/ (p0^1.5) ); E= (-g0/(2*FARADAYS)+s0/(2*FARADAYS)*(T_st-T0)+R_G_U*T_st/(2*FARADAYS).*log(p_H2./p0) + R_G_U*T_st/(4*FARADAYS).*log(p_O2./p0)) ; XX=E-ydata(1:length(E))/1000;%-i_current_density.*tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)) ; %XXAV=E-Average_voltage/1000;%-i_current_density.*tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)) ; XXMED=E- Median_voltage(1:length(E))/1000; XXMIN=E-Min_voltage(1:length(E))/1000; XXMAX=E-Max_voltage(1:length(E))/1000; %[ K1=1/alpha_c; K2=i0 ; K3=pressure_coeff; K4 ; R_membrane;] %[ K1; K2=i0 ; K3=pressure_coeff; K4 ; R_membrane;] % same as (b) model shown in validation section Anotw=[R_G_U*T_st/FARADAYS.*(log(i_app+iloss) +Ec/R_G_U*(1./T_st-1/T0) ) ,-R_G_U*T_st/FARADAYS ,-R_G_U*T_st/FARADAYS.*log(p_O2/p0), tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)).* (i_app+iloss)];%,i_current_density % AnotwMAX=Anotw; AnotwMIN=Anotw; % % %,log(p_O2/p0) % %t_w=.2%0.016; % AN_GDL_surf_blockage = m_l_an*(10000*1000 /(FC_NUMBEROFCELLS* WATERDENSITY * t_w * FC_ACTIVEAREA ) ); % [% blocked] % % AN_GDL_surf_blockage(AN_GDL_surf_blockage >= 1) = 1; % AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - AN_GDL_surf_blockage); %%%IGNORE TWL % %Calculate the apparent current density using the stack current. % i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] % % AnotwMIN=[i_current_density.*tm./lambda_membr .*exp(-b2*(1/303-1./T_st)), i_current_density,R_G_U*T_st/FARADAYS.*(log(i_app+iloss) ),-R_G_U*T_st/FARADAYS.*log(p_O2/p0),-R_G_U*T_st/FARADAYS,ones(size(T_st))]; % % % %t_w=.04%0.02; % %t_w=0.00001; % % AN_GDL_surf_blockage = m_l_an*(10000*1000 /(FC_NUMBEROFCELLS* WATERDENSITY * t_w * FC_ACTIVEAREA ) ); % [% blocked] % % AN_GDL_surf_blockage(AN_GDL_surf_blockage >= 1) = 1; % % AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - AN_GDL_surf_blockage); %%%IGNORE TWL % % % %Calculate the apparent current density using the stack current. % i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] % AnotwMAX=[i_current_density.*tm./lambda_membr .*exp(-b2*(1/303-1./T_st)), i_current_density,R_G_U*T_st/FARADAYS.*(log(i_app+iloss) ),-R_G_U*T_st/FARADAYS.*log(p_O2/p0),-R_G_U*T_st/FARADAYS,ones(size(T_st))]; % %index=find(y(:,7)>.1); params=Anotw(index,:)\XX(index); paramsMED=Anotw(index,:)\XXMED(index); paramsMIN=AnotwMIN(index,:)\XXMIN(index); paramsMAX=AnotwMAX(index,:)\XXMAX(index); % % figure(188) % %plot(t,ydata,'--k',t,1000*(E-i_current_density.*tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st))-Anotw*params),'r:x',time,Min_voltage,'g-.',... % plot(t,ydata(1:length(t)),'--k',t,1000*(E-Anotw*params),'r:x',time,Min_voltage,'g-.',... % time,data( : , 17 ),time,data( : , 2 ),time,data( : , 3 ),time,data( : , 18 ),time,data( : , 19 )) % legend('Average Measured Voltage','Predicted Voltage','Min Voltage') % title('Voltage Model') % ylabel('voltage (mV)') % xlabel('time (s)') figure(6);clf subplot(411); plot(t/60,data(1:LL,2:25),'LineWidth',0.25);hold on % plot(t/60,data(1:LL,3),'LineWidth',0.25,'Color',[0.81 0.81 0.81]);hold on % plot(t/60,data(1:LL,4),'LineWidth',0.25,'Color',[0.89 0.82 0.7]);hold on % plot(t/60,data(1:LL,5),'LineWidth',0.25,'Color',[0.7 0.78 0.87]);hold on % plot(t/60,data(1:LL,6),'LineWidth',0.25,'Color',[0.84 0.77 0.77]);hold on % plot(t/60,data(1:LL,7),'LineWidth',0.25,'Color',[0.74 0.81 0.75]);hold on % plot(t/60,data(1:LL,8),'LineWidth',0.25,'Color',[0.75 0.8 0.8]);hold on % plot(t/60,data(1:LL,9),'LineWidth',0.25,'Color',[0.77 0.82 0.82]);hold on % plot(t/60,data(1:LL,10),'LineWidth',0.25,'Color',[0.77 0.85 0.85]);hold on % plot(t/60,data(1:LL,11),'LineWidth',0.25,'Color',[0.8 0.77 0.85]);hold on % plot(t/60,data(1:LL,12),'LineWidth',0.25,'Color',[0.8 0.75 0.81]);hold on % plot(t/60,data(1:LL,13),'LineWidth',0.25,'Color',[0.85 0.75 0.8]);hold on % plot(t/60,data(1:LL,14),'LineWidth',0.25,'Color',[0.85 0.75 0.85]);hold on % plot(t/60,data(1:LL,15),'LineWidth',0.25,'Color',[0.89 0.87 0.85]);hold on % plot(t/60,data(1:LL,16),'LineWidth',0.25,'Color',[0.89 0.87 0.81]);hold on % plot(t/60,data(1:LL,17),'LineWidth',0.25,'Color',[0.89 0.87 0.75]);hold on % plot(t/60,data(1:LL,18),'LineWidth',0.25,'Color',[0.87 0.8 0.75]);hold on % plot(t/60,data(1:LL,19),'LineWidth',0.25,'Color',[0.75 0.75 0.75]);hold on % plot(t/60,data(1:LL,20),'LineWidth',0.25,'Color',[0.7 0.75 0.7]);hold on % plot(t/60,data(1:LL,21),'LineWidth',0.25,'Color',[0.73 0.8 0.75]);hold on % plot(t/60,data(1:LL,22),'LineWidth',0.25,'Color',[0.7 0.8 0.8]);hold on % plot(t/60,data(1:LL,23),'LineWidth',0.25,'Color',[0.77 0.83 0.77]);hold on % plot(t/60,data(1:LL,24),'LineWidth',0.25,'Color',[0.8 0.7 0.7]);hold on % plot(t/60,data(1:LL,25),'LineWidth',0.25,'Color',[0.9 0.7 0.7]);hold on plot(t/60,ydata(1:length(t)),'b','LineWidth',3);hold on plot(t/60,1000*(E-Anotw*params),'r:','LineWidth',3); %legend('Measured Avg','Estimated Avg') axis([180 200 640 720]) ylabel('Cell Voltages (mV)');xlabel('Time (min)') subplot(412);plot(t/60,i_current_density,'b:','LineWidth',2);hold on; plot(t/60,i_app,'k') ylabel('i (A/cm^2)');xlabel('Time (min)') legend('i','i_{app}') axis([180 200 0.25 0.35]) %subplot(413);[haxes,hline1,hline2]=plotyy(tr/60,ur(:,1)/100000,tr/60,ur(:,4)/100000,'plot','plot') %xlabel('Time (min)'); %axes(haxes(1));ylabel('p_{ca,in} (Bar)');set(hline1,'LineStyle','-','Color','g','LineWidth',2) %axes(haxes(2));ylabel('p_{an,in} (Bar)');set(hline2,'LineStyle',':','Color','r','LineWidth',2) subplot(413);plot(time/60,ur(:,1)/100000,'g:','LineWidth',2);hold on; plot(time/60,ur(:,4)/100000,'r','LineWidth',1) xlabel('Time (min)');ylabel('p (Bar)');legend('Ca In','An In') axis([180 200 1.0 1.25]) subplot(414);plot(time/60,ur(:,15)-273.15,'LineWidth',1) xlabel('Time (min)');ylabel('T (C)') axis([180 200 50 62]) MeanError=mean(abs(ydata(1:length(t))-1000*(E-Anotw*params))./ydata(1:length(t))) MaxError=max(abs(ydata(1:length(t))-1000*(E-Anotw*params))./ydata(1:length(t))) StdError=std(ydata(1:length(t))-1000*(E-Anotw*params)) figure(190);clf subplot(321); plot(t/60,data(1:LL,2:25),'LineWidth',0.25);hold on % plot(t/60,data(1:LL,3),'LineWidth',0.25,'Color',[0.81 0.81 0.81]);hold on % plot(t/60,data(1:LL,4),'LineWidth',0.25,'Color',[0.89 0.82 0.7]);hold on % plot(t/60,data(1:LL,5),'LineWidth',0.25,'Color',[0.7 0.78 0.87]);hold on % plot(t/60,data(1:LL,6),'LineWidth',0.25,'Color',[0.84 0.77 0.77]);hold on % plot(t/60,data(1:LL,7),'LineWidth',0.25,'Color',[0.74 0.81 0.75]);hold on % plot(t/60,data(1:LL,8),'LineWidth',0.25,'Color',[0.75 0.8 0.8]);hold on % plot(t/60,data(1:LL,9),'LineWidth',0.25,'Color',[0.77 0.82 0.82]);hold on % plot(t/60,data(1:LL,10),'LineWidth',0.25,'Color',[0.77 0.85 0.85]);hold on % plot(t/60,data(1:LL,11),'LineWidth',0.25,'Color',[0.8 0.77 0.85]);hold on % plot(t/60,data(1:LL,12),'LineWidth',0.25,'Color',[0.8 0.75 0.81]);hold on % plot(t/60,data(1:LL,13),'LineWidth',0.25,'Color',[0.85 0.75 0.8]);hold on % plot(t/60,data(1:LL,14),'LineWidth',0.25,'Color',[0.85 0.75 0.85]);hold on % plot(t/60,data(1:LL,15),'LineWidth',0.25,'Color',[0.89 0.87 0.85]);hold on % plot(t/60,data(1:LL,16),'LineWidth',0.25,'Color',[0.89 0.87 0.81]);hold on % plot(t/60,data(1:LL,17),'LineWidth',0.25,'Color',[0.89 0.87 0.75]);hold on % plot(t/60,data(1:LL,18),'LineWidth',0.25,'Color',[0.87 0.8 0.75]);hold on % plot(t/60,data(1:LL,19),'LineWidth',0.25,'Color',[0.75 0.75 0.75]);hold on % plot(t/60,data(1:LL,20),'LineWidth',0.25,'Color',[0.7 0.75 0.7]);hold on % plot(t/60,data(1:LL,21),'LineWidth',0.25,'Color',[0.73 0.8 0.75]);hold on % plot(t/60,data(1:LL,22),'LineWidth',0.25,'Color',[0.7 0.8 0.8]);hold on % plot(t/60,data(1:LL,23),'LineWidth',0.25,'Color',[0.77 0.83 0.77]);hold on % plot(t/60,data(1:LL,24),'LineWidth',0.25,'Color',[0.8 0.7 0.7]);hold on % plot(t/60,data(1:LL,25),'LineWidth',0.25,'Color',[0.9 0.7 0.7]);hold on plot(t/60,ydata(1:length(t)),'b','LineWidth',3);hold on plot(t/60,1000*(E-Anotw*params),'r:','LineWidth',3); %legend('Measured Avg','Estimated Avg') axis([183.35 183.6 640 720]) ylabel('Cell Voltages (mV)');xlabel('Time (min)') subplot(323);plot(t/60,i_current_density,'b:','LineWidth',2);hold on; plot(t/60,i_app,'k') ylabel('i (A/cm^2)');xlabel('Time (min)') legend('i','i_{app}') axis([183.35 183.6 0.25 0.35]) %subplot(325);plot(tr/60,ur(:,1)/100000,'g',tr/60,ur(:,4)/100000,'r:') %xlabel('Time (min)');ylabel('p (Bar)');legend('Ca In','An In') %axis([183.35 183.45 1 1.4]) subplot(325);plot(tout/60,p_o2(:,4),'c');hold on; plot(tout/60,p_o2(:,3),'g-.');hold on; plot(tout/60,p_o2(:,2),'r:');hold on; plot(tout/60,p_o2(:,1),'b--') %axis([0 944 0 .13]) xlabel('Time (min)');ylabel('p_{O_2} (bar)');legend('ch','3','2','mb'); axis([183.35 183.6 0.05 0.15]) subplot(322); plot(t/60,data(1:LL,2:25),'LineWidth',0.25);hold on%,'Color',[0.8 0.9 0.9]);hold on % plot(t/60,data(1:LL,3),'LineWidth',0.25,'Color',[0.81 0.81 0.81]);hold on % plot(t/60,data(1:LL,4),'LineWidth',0.25,'Color',[0.89 0.82 0.7]);hold on % plot(t/60,data(1:LL,5),'LineWidth',0.25,'Color',[0.7 0.78 0.87]);hold on % plot(t/60,data(1:LL,6),'LineWidth',0.25,'Color',[0.84 0.77 0.77]);hold on % plot(t/60,data(1:LL,7),'LineWidth',0.25,'Color',[0.74 0.81 0.75]);hold on % plot(t/60,data(1:LL,8),'LineWidth',0.25,'Color',[0.75 0.8 0.8]);hold on % plot(t/60,data(1:LL,9),'LineWidth',0.25,'Color',[0.77 0.82 0.82]);hold on % plot(t/60,data(1:LL,10),'LineWidth',0.25,'Color',[0.77 0.85 0.85]);hold on % plot(t/60,data(1:LL,11),'LineWidth',0.25,'Color',[0.8 0.77 0.85]);hold on % plot(t/60,data(1:LL,12),'LineWidth',0.25,'Color',[0.8 0.75 0.81]);hold on % plot(t/60,data(1:LL,13),'LineWidth',0.25,'Color',[0.85 0.75 0.8]);hold on % plot(t/60,data(1:LL,14),'LineWidth',0.25,'Color',[0.85 0.75 0.85]);hold on % plot(t/60,data(1:LL,15),'LineWidth',0.25,'Color',[0.89 0.87 0.85]);hold on % plot(t/60,data(1:LL,16),'LineWidth',0.25,'Color',[0.89 0.87 0.81]);hold on % plot(t/60,data(1:LL,17),'LineWidth',0.25,'Color',[0.89 0.87 0.75]);hold on % plot(t/60,data(1:LL,18),'LineWidth',0.25,'Color',[0.87 0.8 0.75]);hold on % plot(t/60,data(1:LL,19),'LineWidth',0.25,'Color',[0.75 0.75 0.75]);hold on % plot(t/60,data(1:LL,20),'LineWidth',0.25,'Color',[0.7 0.75 0.7]);hold on % plot(t/60,data(1:LL,21),'LineWidth',0.25,'Color',[0.73 0.8 0.75]);hold on % plot(t/60,data(1:LL,22),'LineWidth',0.25,'Color',[0.7 0.8 0.8]);hold on % plot(t/60,data(1:LL,23),'LineWidth',0.25,'Color',[0.77 0.83 0.77]);hold on % plot(t/60,data(1:LL,24),'LineWidth',0.25,'Color',[0.8 0.7 0.7]);hold on % plot(t/60,data(1:LL,25),'LineWidth',0.25,'Color',[0.9 0.7 0.7]);hold on plot(t/60,ydata(1:length(t)),'b','LineWidth',3); plot(t/60,1000*(E-Anotw*params),'r:','LineWidth',3); %legend('Measured Avg','Estimated Avg') axis([198.22 198.4 650 720]) ylabel('Cell Voltages (mV)');xlabel('Time (min)') subplot(324);plot(t/60,i_current_density,'b:','LineWidth',2);hold on plot(t/60,i_app,'k') ylabel('i (A/cm^2)');xlabel('Time (min)') legend('i','i_{app}') axis([198.22 198.4 0.25 0.35]) %subplot(326);plot(tr/60,ur(:,1)/100000,'g',tr/60,ur(:,4)/100000,'r:') %xlabel('Time (min)');ylabel('p (Bar)');legend('Ca In','An In') %axis([198.2 198.3 1 1.4]) subplot(326);plot(tout/60,p_o2(:,4),'c','LineWidth',2);hold on; plot(tout/60,p_o2(:,3),'g-.','LineWidth',2);hold on; plot(tout/60,p_o2(:,2),'r:','LineWidth',2);hold on; plot(tout/60,p_o2(:,1),'b--','LineWidth',2) xlabel('Time (min)');ylabel('p_{O_2} (bar)');legend('ch','3','2','mb'); axis([198.22 198.4 0.05 0.15]) % plots for statisitcal based models %time,Min_voltage,t,1000*(E-AnotwMIN*paramsMIN),'g--',... % time,Max_voltage,t,1000*(E-AnotwMAX*paramsMAX),'--k',time,Median_voltage, t,1000*(E-Anotw*paramsMED),':',... % time(index),605*ones(size(index)),'kd') % time,data( : , 17 ),time,data( : , 2 ),time,data( : , 3 ),time,data( : , 18 ),time,data( : , 19 )) figure(192);clf subplot(421); plot(t/60,data(1:LL,2:25),'LineWidth',0.25);hold on plot(t/60,ydata(1:length(t)),'b','LineWidth',3);hold on plot(t/60,1000*(E-Anotw*params),'r:','LineWidth',3); %legend('Measured Avg','Estimated Avg') axis([183.2 184 640 720]) ylabel('Cell Voltages (mV)');xlabel('Time (min)') subplot(322);plot(t/60,i_current_density,'b:','LineWidth',2);hold on; plot(t/60,i_app,'k') ylabel('i (A/cm^2)');xlabel('Time (min)') legend('i','i_{app}') axis([183.2 184 0.25 0.35]) %subplot(325);plot(tr/60,ur(:,1)/100000,'g',tr/60,ur(:,4)/100000,'r:') %xlabel('Time (min)');ylabel('p (Bar)');legend('Ca In','An In') %axis([183.35 183.45 1 1.4]) subplot(423);plot(tout/60,p_o2(:,4),'c');hold on; plot(tout/60,p_o2(:,3),'g-.');hold on; plot(tout/60,p_o2(:,2),'r:');hold on; plot(tout/60,p_o2(:,1),'b--') %axis([0 944 0 .13]) xlabel('Time (min)');ylabel('p_{O_2} (bar)');legend('4','3','2','1'); axis([183.2 184 0.05 0.15]) subplot(424);plot(tout/60,p_h2(:,4),'c','LineWidth',1);hold on; plot(tout/60,p_h2(:,3),'g-.','LineWidth',1);hold on; plot(tout/60,p_h2(:,2),'r:','LineWidth',1);hold on; plot(tout/60,p_h2(:,1),'b--','LineWidth',1) xlabel('Time (min)');ylabel('p_{H_2} (bar)');legend('4','3','2','1'); axis([183.2 184 0.95 1.15]) subplot(425); %figure plot(tout/60,CATHODE_Pressure(:,2)/100000,'g-.','LineWidth',2);hold on; plot(tout/60,CATHODE_3(:,9).*R_G_U.*CATHODE_1(:,3)/100000,'g-.','LineWidth',2);hold on; plot(tout/60,CATHODE_3(:,9).*R_G_U.*CATHODE_2(:,3)/100000,'r:','LineWidth',2);hold on; plot(tout/60,CATHODE_3(:,9).*R_G_U.*CATHODE_3(:,3)/100000,'b--','LineWidth',2) xlabel('Time (min)');ylabel('p_{v,ca} (bar)');legend('4','3','2','1'); axis([183.2 184 0.13 0.16]) subplot(426); plot(tout/60,p_v_an(:,4),'c','LineWidth',2);hold on; plot(tout/60,p_v_an(:,3),'g-.','LineWidth',2);hold on; plot(tout/60,p_v_an(:,2),'r:','LineWidth',2);hold on; plot(tout/60,p_v_an(:,1),'b--','LineWidth',2) xlabel('Time (min)');ylabel('p_{v,an} (bar)');legend('4','3','2','1'); axis([183.2 184 0.125 0.145]) subplot(427); plot(tout/60,CATHODE_1(:,7),'g-.','LineWidth',2);hold on; plot(tout/60,CATHODE_2(:,7),'r:','LineWidth',2);hold on; plot(tout/60,CATHODE_3(:,7),'b--','LineWidth',2) xlabel('Time (min)');ylabel('s_{ca}');legend('3','2','1'); axis([183.2 184 0.13 0.16]) subplot(428); plot(tout/60,ANODE_1(:,5),'g-.','LineWidth',2);hold on; plot(tout/60,ANODE_2(:,5),'r:','LineWidth',2);hold on; plot(tout/60,ANODE_3(:,5),'b--','LineWidth',2) xlabel('Time (min)');ylabel('s_{an}');legend('3','2','1'); axis([183.2 184 0.1 0.14]) figure(193);clf subplot(221); colorsp=get(gca,'ColorOrder') colorspn=rgb2hsv(colorsp); colorspn(:,2)=colorspn(:,2)*.3; %xx=rand([1,24]) colorspn=[linspace(0,1,24).' , [0.6787 0.7577 0.7431 0.3922 0.6555 0.1712 0.7060 0.0318 0.2769 0.0462 0.0971 0.8235 0.6948 0.3171 0.9502 0.0344 0.4387 0.3816 0.7655 0.7952 0.1869 0.4898 0.4456 0.6463] , 0.8*ones(24,1)]% *.4 colorspn=hsv2rgb(colorspn); set(gca,'ColorOrder',colorspn,'NextPlot','add')%'replacechildren') plot(t/60,data(1:LL,2:25),'LineWidth',0.25);hold on plot(t/60,ydata(1:length(t)),'b','LineWidth',2.5);hold on plot(t/60,1000*(E-Anotw*params),'r-.','LineWidth',2.5); %legend('Measured Avg','Estimated Avg') axis([183.2 184 640 720]) ylabel('Cell Voltages (mV)');xlabel('Time (min)') box on subplot(222);plot(t/60,i_current_density,'b:','LineWidth',2);hold on; plot(t/60,i_app,'k') ylabel('i (A/cm^2)');xlabel('Time (min)') legend('i','i_{app}') axis([183.2 184 0.25 0.35]) %subplot(325);plot(tr/60,ur(:,1)/100000,'g',tr/60,ur(:,4)/100000,'r:') %xlabel('Time (min)');ylabel('p (Bar)');legend('Ca In','An In') %axis([183.35 183.45 1 1.4]) subplot(223);plot(tout/60,p_o2(:,4),'c','LineWidth',1.5);hold on; plot(tout/60,p_o2(:,3),'g-.','LineWidth',1.5);hold on; plot(tout/60,p_o2(:,2),'r:','LineWidth',1.5);hold on; plot(tout/60,p_o2(:,1),'b--','LineWidth',1.5) %axis([0 944 0 .13]) xlabel('Time (min)');ylabel('p_{O_2} (bar)');legend('4','3','2','1'); axis([183.2 184 0.05 0.15]) subplot(224);plot(tout/60,p_h2(:,4),'c','LineWidth',1.5);hold on; plot(tout/60,p_h2(:,3),'g-.','LineWidth',1.5);hold on; plot(tout/60,p_h2(:,2),'r:','LineWidth',1.5);hold on; plot(tout/60,p_h2(:,1),'b--','LineWidth',1.5) xlabel('Time (min)');ylabel('p_{H_2} (bar)');legend('4','3','2','1'); axis([183.2 184 0.95 1.15]) %print -depsc -r300 figure193 figure(194);clf subplot(321); %figure plot(tout/60,CATHODE_Pressure(:,2)/100000,'c','LineWidth',1.5);hold on; plot(tout/60,CATHODE_3(:,9).*R_G_U.*CATHODE_1(:,3)/100000,'g-.','LineWidth',1.5);hold on; plot(tout/60,CATHODE_3(:,9).*R_G_U.*CATHODE_2(:,3)/100000,'r:','LineWidth',1.5);hold on; plot(tout/60,CATHODE_3(:,9).*R_G_U.*CATHODE_3(:,3)/100000,'b--','LineWidth',1.5) xlabel('Time (min)');ylabel('p_{v,ca} (bar)');legend('4','3','2','1'); axis([183.2 184 0.13 0.16]) subplot(322); plot(tout/60,p_v_an(:,4),'c','LineWidth',1.5);hold on; plot(tout/60,p_v_an(:,3),'g-.','LineWidth',1.5);hold on; plot(tout/60,p_v_an(:,2),'r:','LineWidth',1.5);hold on; plot(tout/60,p_v_an(:,1),'b--','LineWidth',1.5) xlabel('Time (min)');ylabel('p_{v,an} (bar)');legend('4','3','2','1'); axis([183.2 184 0.125 0.145]) subplot(323); plot(tout/60,CATHODE_1(:,7),'g-.','LineWidth',1.5);hold on; plot(tout/60,CATHODE_2(:,7),'r:','LineWidth',1.5);hold on; plot(tout/60,CATHODE_3(:,7),'b--','LineWidth',1.5) xlabel('Time (min)');ylabel('s_{ca}');legend('3','2','1'); axis([183.2 184 0.13 0.16]) subplot(324); plot(tout/60,ANODE_1(:,5),'g-.','LineWidth',1.5);hold on; plot(tout/60,ANODE_2(:,5),'r:','LineWidth',1.5);hold on; plot(tout/60,ANODE_3(:,5),'b--','LineWidth',1.5) xlabel('Time (min)');ylabel('s_{an}');legend('3','2','1'); axis([183.2 184 0.1 0.14]) subplot(325); plot(time/60,ur(:,16),'LineWidth',1.5); xlabel('Time (min)');ylabel('Purge State'); axis([183.2 184 -0.1 1.1]) subplot(326); plot(tout/60,1000*y(:,8),'LineWidth',1.5) xlabel('Time (min)');ylabel('m_{lw,an} (mg)'); axis([183.2 184 0 5]) % print -depsc -r300 figure194 % print -depsc -r100 figure194a1 K1=params(1) K2=exp(params(2)/params(1)) K3=params(3)/params(1) K4=params(4) %Rohm=params(5) K1=paramsMED(1) K2=exp(paramsMED(2)/paramsMED(1)) K3=paramsMED(3)/paramsMED(1) K4=paramsMED(4) %Rohm=paramsMED(5) K1=paramsMIN(1) K2=exp(paramsMIN(2)/paramsMIN(1)) K3=paramsMIN(3)/paramsMIN(1) K4=paramsMIN(4) %Rohm=paramsMIN(5) % % load Fit_Model % % % figure(190) % plot(t,ydata(1:length(t)),t,1000*(E-Anotw*params),'r-.',time,Min_voltage,t,1000*(E-AnotwMIN*paramsMIN),'g--',... % time,Max_voltage,t,1000*(E-AnotwMAX*paramsMAX),'--k',time,Median_voltage, t,1000*(E-Anotw*paramsMED),':') % % time,data( : , 17 ),time,data( : , 2 ),time,data( : , 3 ),time,data( : , 18 ),time,data( : , 19 )) % legend('Average Measured Voltage','Predicted Voltage (AVG)','Min Voltage','Predicted Voltage (MIN)','Max Voltage','Predicted Voltage (MAX)','Median Voltage','Predicted Voltage (MED)') % title('Voltage Model Validation') % ylabel('voltage (mV)') % xlabel('time (s)') t_wa =.13146%/1.2%/2; t_wb =.13146%*1.2; t_wc =.13146%/1.3; AN_GDL_surf_blockage = m_l_an*(10000*1000 /(FC_NUMBEROFCELLS* WATERDENSITY * t_wa * FC_ACTIVEAREA ) ); % [% blocked] AN_GDL_surf_blockage(AN_GDL_surf_blockage >= 1) = 1; AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - 2*AN_GDL_surf_blockage); %%%IGNORE TWL %Calculate the apparent current density using the stack current. i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] Anotwa=[R_G_U*T_st/FARADAYS.*(log(i_app+iloss) +Ec/R_G_U*(1./T_st-1/T0) ) ,-R_G_U*T_st/FARADAYS ,-R_G_U*T_st/FARADAYS.*log(p_O2/p0), tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)).* (i_current_density+iloss)]; paramsa=Anotwa(index,:)\XX(index); AN_GDL_surf_blockage = m_l_an*(10000*1000 /(FC_NUMBEROFCELLS* WATERDENSITY * t_wb * FC_ACTIVEAREA ) ); % [% blocked] AN_GDL_surf_blockage(AN_GDL_surf_blockage >= 1) = 1; AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - 2*AN_GDL_surf_blockage); %%%IGNORE TWL %Calculate the apparent current density using the stack current. i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] Anotwb=[R_G_U*T_st/FARADAYS.*(log(i_app+iloss) +Ec/R_G_U*(1./T_st-1/T0) ) ,-R_G_U*T_st/FARADAYS ,-R_G_U*T_st/FARADAYS.*log(p_O2/p0), tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)).* (i_app+iloss)]; paramsb=Anotwb(index,:)\XX(index); AN_GDL_surf_blockage = m_l_an*(10000*1000 /(FC_NUMBEROFCELLS* WATERDENSITY * t_wc * FC_ACTIVEAREA ) ); % [% blocked] AN_GDL_surf_blockage(AN_GDL_surf_blockage >= 1) = 1; AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - 2*AN_GDL_surf_blockage); %%%IGNORE TWL %Calculate the apparent current density using the stack current. i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] Anotwc=[R_G_U*T_st/FARADAYS.*(log(i_app+iloss) +Ec/R_G_U*(1./T_st-1/T0) ) ,-R_G_U*T_st/FARADAYS ,-R_G_U*T_st/FARADAYS.*log(p_O2/p0), tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)).* (i_app+iloss),i_current_density]; paramsc=Anotwc(index,:)\XX(index); % K1=paramsa(1) % K2=exp(paramsa(2)/paramsa(1)) % K3=paramsa(3)/paramsa(1) % K4=paramsa(4) %%%%%%%%%%%%%%%%%%%%%% % voltage parameters back calculated for voltage model (b) K1=paramsb(1) K2=exp(paramsb(2)/paramsb(1)) K3=paramsb(3)/paramsb(1) K4=paramsb(4) % K1=paramsc(1) % K2=exp(paramsc(2)/paramsc(1)) % K3=paramsc(3)/paramsc(1) % K4=paramsc(4) % Rohm=paramsc(5) % figure(192) % plot(t,ydata(1:length(t)),t,1000*(E-Anotwa*paramsa),'r-.',t,1000*(E-Anotwb*paramsb),'g--',t,1000*(E-Anotwc*paramsc),'k:',... % time(index),605*ones(size(index)),'kd') % legend('Average Measured Voltage','Predicted Voltage (AVGa)','Predicted Voltage (AVGb)','Predicted Voltage (AVGc)') % title('Voltage Calibration Data (comparison)') % ylabel('voltage (mV)') % xlabel('time (s)') % % figure(193) % plot(t,ydata(1:length(t)),t,1000*(E-Anotwb*paramsb),'g--',time,Min_voltage,time,Max_voltage,... % time(index),605*ones(size(index)),'kd') % legend('Average Measured Voltage','Predicted Voltage','representative min','representative max') % title('Voltage Calibration Data') % ylabel('voltage (mV)') % xlabel('time (s)') %return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Validation Data %load Dec17_validation_alpha_11_5 load Dec17_validation_alpha_11_5_all mscaledata_a; %load Tuning_data_Dec21_1109 % xdata_valid=y(:,[2 3 4 6 7 8 ]); % time_valid=t; % data_valid=data; clear AveV; for i=1:length(data(:,1)) AveV(i)=sum(data(i,2:25))/24; end T_st=xdata_valid(:,1); % [K] c_O2=xdata_valid(:,2); % [pa] c_H2=xdata_valid(:,3); % [pa] p_O2=c_O2.*T_st*R_G_U; % calculate pressure from the ideal gas law p_H2=c_H2.*T_st*R_G_U; lambda_membr=xdata_valid(:,4); I_st=xdata_valid(:,5); % [A] m_l_an=xdata_valid(:,6); % [kg] i_current_density=I_st./FC_ACTIVEAREA; %[A/cm^2] AN_GDL_surf_blockage = m_l_an*(10000*1000 /(FC_NUMBEROFCELLS* WATERDENSITY * t_wa * FC_ACTIVEAREA ) ); % [% blocked] AN_GDL_surf_blockage(AN_GDL_surf_blockage >= 1) = 1; AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - 2*AN_GDL_surf_blockage); %%%IGNORE TWL %Calculate the apparent current density using the stack current. i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] Anotwa=[R_G_U*T_st/FARADAYS.*(log(i_app+iloss) +Ec/R_G_U*(1./T_st-1/T0) ) ,-R_G_U*T_st/FARADAYS ,-R_G_U*T_st/FARADAYS.*log(p_O2/p0), tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)).* (i_current_density+iloss)]; %paramsa=Anotwa(index,:)\XX(index); AN_GDL_surf_blockage = m_l_an*(10000*1000 /(FC_NUMBEROFCELLS* WATERDENSITY * t_wb * FC_ACTIVEAREA ) ); % [% blocked] AN_GDL_surf_blockage(AN_GDL_surf_blockage >= 1) = 1; AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - 2*AN_GDL_surf_blockage); %%%IGNORE TWL %Calculate the apparent current density using the stack current. i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] Anotwb=[R_G_U*T_st/FARADAYS.*(log(i_app+iloss) +Ec/R_G_U*(1./T_st-1/T0) ) ,-R_G_U*T_st/FARADAYS ,-R_G_U*T_st/FARADAYS.*log(p_O2/p0), tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)).* (i_app+iloss)]; %paramsb=Anotwb(index,:)\XX(index); AN_GDL_surf_blockage = m_l_an*(10000*1000 /(FC_NUMBEROFCELLS* WATERDENSITY * t_wc * FC_ACTIVEAREA ) ); % [% blocked] AN_GDL_surf_blockage(AN_GDL_surf_blockage >= 1) = 1; AN_fc_activearea_eff = FC_ACTIVEAREA * (1 - 2*AN_GDL_surf_blockage); %%%IGNORE TWL %Calculate the apparent current density using the stack current. i_app=I_st./AN_fc_activearea_eff; %[A/cm^2] Anotwc=[R_G_U*T_st/FARADAYS.*(log(i_app+iloss) +Ec/R_G_U*(1./T_st-1/T0) ) ,-R_G_U*T_st/FARADAYS ,-R_G_U*T_st/FARADAYS.*log(p_O2/p0), tm./(b11.*lambda_membr-b12) .*exp(-b2*(1/303-1./T_st)).* (i_app+iloss),i_current_density]; %paramsc=Anotwc(index,:)\XX(index); E= (-g0/(2*FARADAYS)+s0/(2*FARADAYS)*(T_st-T0)+R_G_U*T_st/(2*FARADAYS).*log(p_H2./p0) + R_G_U*T_st/(4*FARADAYS).*log(p_O2./p0)) ; % figure(194) % plot(tr_valid,ydata_valid(1:length(tr_valid)),tr_valid,1000*(E-Anotwa*paramsa),'r-.',tr_valid,1000*(E-Anotwb*paramsb),'g--',tr_valid,1000*(E-Anotwc*paramsc),'k:') % legend('Average Measured Voltage','Predicted Voltage (AVGa)','Predicted Voltage (AVGb)','Predicted Voltage (AVGc)') % title('Voltage Validation Data ') % ylabel('voltage (mV)') % xlabel('time (s)') % % % % figure(195) % plot(time_valid,data_valid( : , 2:(2+23) ),tr_valid,1000*(E-Anotwb*paramsb),'k--') % %legend('Average Measured Voltage','Predicted Voltage (AVGa)','Predicted Voltage (AVGb)','Predicted Voltage (AVGc)') % title('Voltage Validation Data ') % ylabel('voltage (mV)') % xlabel('time (s)') % % % figure(196); plot(tr_valid,i_current_density,tr_valid,i_app,'k-.') % ylabel('Current Density ( A / cm^2 )') % xlabel('Time (s)') % legend('i','i_{app}') % title('Validation Data') figure(300);clf subplot(411); plot(time_valid/60,data_valid(:,2:25));hold on%,'LineWidth',0.25,'Color',[0.8 0.9 0.9]);hold on % plot(time_valid/60,data_valid(:,3),'LineWidth',0.25,'Color',[0.81 0.81 0.81]);hold on % plot(time_valid/60,data_valid(:,4),'LineWidth',0.25,'Color',[0.89 0.82 0.7]);hold on % plot(time_valid/60,data_valid(:,5),'LineWidth',0.25,'Color',[0.7 0.78 0.87]);hold on % plot(time_valid/60,data_valid(:,6),'LineWidth',0.25,'Color',[0.84 0.77 0.77]);hold on % plot(time_valid/60,data_valid(:,7),'LineWidth',0.25,'Color',[0.74 0.81 0.75]);hold on % plot(time_valid/60,data_valid(:,8),'LineWidth',0.25,'Color',[0.75 0.8 0.8]);hold on % plot(time_valid/60,data_valid(:,9),'LineWidth',0.25,'Color',[0.77 0.82 0.82]);hold on % plot(time_valid/60,data_valid(:,10),'LineWidth',0.25,'Color',[0.77 0.85 0.85]);hold on % plot(time_valid/60,data_valid(:,11),'LineWidth',0.25,'Color',[0.8 0.77 0.85]);hold on % plot(time_valid/60,data_valid(:,12),'LineWidth',0.25,'Color',[0.8 0.75 0.81]);hold on % plot(time_valid/60,data_valid(:,13),'LineWidth',0.25,'Color',[0.85 0.75 0.8]);hold on % plot(time_valid/60,data_valid(:,14),'LineWidth',0.25,'Color',[0.85 0.75 0.85]);hold on % plot(time_valid/60,data_valid(:,15),'LineWidth',0.25,'Color',[0.89 0.87 0.85]);hold on % plot(time_valid/60,data_valid(:,16),'LineWidth',0.25,'Color',[0.89 0.87 0.81]);hold on % plot(time_valid/60,data_valid(:,17),'LineWidth',0.25,'Color',[0.89 0.87 0.75]);hold on % plot(time_valid/60,data_valid(:,18),'LineWidth',0.25,'Color',[0.87 0.8 0.75]);hold on % plot(time_valid/60,data_valid(:,19),'LineWidth',0.25,'Color',[0.75 0.75 0.75]);hold on % plot(time_valid/60,data_valid(:,20),'LineWidth',0.25,'Color',[0.7 0.75 0.7]);hold on % plot(time_valid/60,data_valid(:,21),'LineWidth',0.25,'Color',[0.73 0.8 0.75]);hold on % plot(time_valid/60,data_valid(:,22),'LineWidth',0.25,'Color',[0.7 0.8 0.8]);hold on % plot(time_valid/60,data_valid(:,23),'LineWidth',0.25,'Color',[0.77 0.83 0.77]);hold on % plot(time_valid/60,data_valid(:,24),'LineWidth',0.25,'Color',[0.8 0.7 0.7]);hold on % plot(time_valid/60,data_valid(:,25),'LineWidth',0.25,'Color',[0.9 0.7 0.7]);hold on plot(time_valid/60,AveV,'b','LineWidth',3);hold on plot(time_valid/60,1000*(E-Anotwb*paramsb),'r:','LineWidth',3); %legend('Measured Avg','Estimated Avg') %axis([191 208 720 780]) ylabel('Cell Voltages (mV)');xlabel('Time (min)') subplot(412);plot(time_valid/60,i_current_density,'b:','LineWidth',2);hold on plot(time_valid/60,i_app,'k') ylabel('Current Density (A/cm^2)');xlabel('Time (min)') legend('i','i_{app}') %axis([191 208 0.08 0.17]) subplot(413);plot(time_valid/60,ur(:,1)/100000,'g:');hold on; plot(time_valid/60,ur(:,4)/100000,'r') xlabel('Time (min)');ylabel('p (Bar)');legend('Ca In','An In') %axis([191 208 1 1.25]) subplot(414);plot(time_valid/60,T_st-273.15) xlabel('Time (min)');ylabel('T (C)') %axis([191 208 48 52]) MeanError=mean(abs(AveV'-1000*(E-Anotwb*paramsb))./AveV') MaxError=max(abs(AveV'-1000*(E-Anotwb*paramsb))./AveV') StdError=std(AveV'-1000*(E-Anotwb*paramsb)) figure(9);clf subplot(411); plot(time_valid/60,data_valid(:,2:25));hold on%,'LineWidth',0.25,'Color',[0.8 0.9 0.9]);hold on % plot(time_valid/60,data_valid(:,3),'LineWidth',0.25,'Color',[0.81 0.81 0.81]);hold on % plot(time_valid/60,data_valid(:,4),'LineWidth',0.25,'Color',[0.89 0.82 0.7]);hold on % plot(time_valid/60,data_valid(:,5),'LineWidth',0.25,'Color',[0.7 0.78 0.87]);hold on % plot(time_valid/60,data_valid(:,6),'LineWidth',0.25,'Color',[0.84 0.77 0.77]);hold on % plot(time_valid/60,data_valid(:,7),'LineWidth',0.25,'Color',[0.74 0.81 0.75]);hold on % plot(time_valid/60,data_valid(:,8),'LineWidth',0.25,'Color',[0.75 0.8 0.8]);hold on % plot(time_valid/60,data_valid(:,9),'LineWidth',0.25,'Color',[0.77 0.82 0.82]);hold on % plot(time_valid/60,data_valid(:,10),'LineWidth',0.25,'Color',[0.77 0.85 0.85]);hold on % plot(time_valid/60,data_valid(:,11),'LineWidth',0.25,'Color',[0.8 0.77 0.85]);hold on % plot(time_valid/60,data_valid(:,12),'LineWidth',0.25,'Color',[0.8 0.75 0.81]);hold on % plot(time_valid/60,data_valid(:,13),'LineWidth',0.25,'Color',[0.85 0.75 0.8]);hold on % plot(time_valid/60,data_valid(:,14),'LineWidth',0.25,'Color',[0.85 0.75 0.85]);hold on % plot(time_valid/60,data_valid(:,15),'LineWidth',0.25,'Color',[0.89 0.87 0.85]);hold on % plot(time_valid/60,data_valid(:,16),'LineWidth',0.25,'Color',[0.89 0.87 0.81]);hold on % plot(time_valid/60,data_valid(:,17),'LineWidth',0.25,'Color',[0.89 0.87 0.75]);hold on % plot(time_valid/60,data_valid(:,18),'LineWidth',0.25,'Color',[0.87 0.8 0.75]);hold on % plot(time_valid/60,data_valid(:,19),'LineWidth',0.25,'Color',[0.75 0.75 0.75]);hold on % plot(time_valid/60,data_valid(:,20),'LineWidth',0.25,'Color',[0.7 0.75 0.7]);hold on % plot(time_valid/60,data_valid(:,21),'LineWidth',0.25,'Color',[0.73 0.8 0.75]);hold on % plot(time_valid/60,data_valid(:,22),'LineWidth',0.25,'Color',[0.7 0.8 0.8]);hold on % plot(time_valid/60,data_valid(:,23),'LineWidth',0.25,'Color',[0.77 0.83 0.77]);hold on % plot(time_valid/60,data_valid(:,24),'LineWidth',0.25,'Color',[0.8 0.7 0.7]);hold on % plot(time_valid/60,data_valid(:,25),'LineWidth',0.25,'Color',[0.9 0.7 0.7]);hold on plot(time_valid/60,AveV,'b','LineWidth',3);hold on plot(time_valid/60,1000*(E-Anotwb*paramsb),'r:','LineWidth',3); %legend('Measured Avg','Estimated Avg') axis([150 220 650 820]) ylabel('Cell Voltages (mV)');xlabel('Time (min)') subplot(412);plot(time_valid/60,i_current_density,'b:','LineWidth',2);hold on plot(time_valid/60,i_app,'k') ylabel('i (A/cm^2)');xlabel('Time (min)') legend('i','i_{app}') axis([150 220 0 0.25]) subplot(413);plot(time/60,ur(:,1)/100000,'g');hold on; plot(time/60,ur(:,4)/100000,'r:') xlabel('Time (min)');ylabel('p (Bar)');legend('Ca In','An In') axis([150 220 1 1.25]) subplot(414);plot(time_valid/60,T_st-273.15) xlabel('Time (min)');ylabel('T (C)') axis([150 220 48 54])