% truck_simulation_UKACC2018.m % % simulation codes for % "Extremum Seeking Control for Truck Drag Reduction" % G Papageorgiou, J Barden, J Whidborne, K Garry % 12th UKACC International Conference on Control % Sheffield, UK, 5th - 7th September 2018 % (c) JF Whidborne, Cranfield University, 12 June 2018 clear all fs=14; % fontsize for plots Tvideo=0.033 % video sample time %% Vehicle Constants rho = 1.22; % air density (kg/m^3) g = 9.81; % gravitational constant (m/s^2) A = 2.0^2; % truck reference area (m^2) Crr = 6.0e-3; % rolling resistance coeeficient Crv = 10.0e-6; % velocity dependent rolling resistance coeeficient m = 20000; % vehicle mass (kg) %% Cruise controller PID coefficients Kp = 1.5 ; % proportional gain Ki = 0.3 ; % integral gain Ka = 1; % antiwindup gain Th2fE = 12000; % throttle to force actuator(400bhp * 746W / 24m/s) %% Aerodynamics (from Barden files - wtDat) [deflector.deltaH,deflector.beta,deflector.CD]=getCDmap; %% sensing and actuating %sn=1e-6; % sensor noise power slewratemax=0.05; % maximum and minimum deflector height ratio slew rate %% ESC controller coefficients % a = 0.01; % ESC perturbation amplitude am = 1; % ESC modulation amplitude omega = 5; % ESC perturbation frequency (rad/s) Kes = -5; % ESC integrator gain (minimizing hence negative) (-1 convergence is slow) Kw = Kes; % ESC antiwindup gain (minimizing hence negative) wh = .1; % ESC high-pass filter coefficient wl = 0.5; % ESC low-pass filter coefficient thetahatmax = max(deflector.deltaH)-a; % ESC antiwindup maximum thetahatmin = min(deflector.deltaH)+a; % ESC antiwindup miximum %% simulate cases casenumber = 2 % Test case number switch casenumber case 1 % zero yaw, wind gradient, constant speed 24m/s izero = find(deflector.beta==0); % get zero beta [CDmin,imin] = min(deflector.CD(:,izero));% get optimal point deltaHmin = deflector.deltaH(imin); %CD0=6*(thetahatmin-1).^2+.54 VTref = 24; % reference speed m/s VT0 = VTref; % initial speed; fM0 = m*g*(Crr+Crv*VT0^2 ); % initial mechanical force % disturbances VTrefchange1 = 0; % no reference speed change 1 m/s tVTrefchange1 = 20; % reference speed change time 1 s VTrefchange2 = 0; % no reference speed change 2 m/s tVTrefchange2 = 50; % reference speed change time 2 s VW0 = 0; % no wind alphachange = 0; % no gradient talphachange = 0; % no gradient change phichange = 0; % no wind heading change tphichange = 0; % no wind heading change % simulate from min deltaH tfinal = 30; % simulation run time (s) thetahat0 = thetahatmin; % theta initial condition CD0 = (deflector.CD(11,izero)); % initial CD fD0 = VT0^2*rho*CD0*A/2; % initial drag force Ki0 = (fD0+fM0)/Th2fE; % initial PI control output sim('Deflector_esc_UKACC2018'); deltaH_minstart = deltaH; CD_minstart = CD; theta_minstart = theta; % simulate from max deltaH thetahat0 = thetahatmax; % theta initial condition CD0 = (deflector.CD(291,izero)); % initial CD fD0 = VT0^2*rho*CD0*A/2; % initial drag force Ki0 = (fD0+fM0)/Th2fE; % initial PI control output sim('Deflector_esc_UKACC2018'); deltaH_maxstart = deltaH; CD_maxstart = CD; theta_maxstart = theta; % plot results figure(11) plot(deltaH_minstart(:,1),deltaH_minstart(:,2),'-',... deltaH_maxstart(:,1),deltaH_maxstart(:,2),'-',... [0,tfinal],[min(deflector.deltaH),min(deflector.deltaH)],'-.',... [0,tfinal],[max(deflector.deltaH),max(deflector.deltaH)],'-.',... [0,tfinal],[deltaHmin,deltaHmin],'--',... 'linewidth',1.2) grid on h99=gca; %h99.XTick=[0:2:14]; grid on ylabel('$\delta_H$','Interpreter' ,'Latex','fontsize',fs) xlabel('$t$ (s)','Interpreter' ,'Latex','fontsize',fs); % g99=legend('$\theta(t)$','$\delta_H^0$ ' ); % set(g99,'Interpreter' ,'Latex','fontsize',12,'Location','southeast') if(0) % print fig = gcf; fig.PaperPositionMode = 'auto'; fig_pos = fig.PaperPosition; fig.PaperSize = [fig_pos(3) fig_pos(4)]; saveas(gcf,'case1_dH', 'pdf') end %% visualize xx=[theta_minstart(:,2);theta_maxstart(:,2)]; yy=[CD_minstart(:,2);CD_maxstart(:,2)]; [x,xisort]=sort(xx);y=yy(xisort); xt=theta_maxstart(:,2); yt=CD_maxstart(:,2); Visualize_esc(x,y,xt,yt); if(0) % print fig = gcf; fig.PaperPositionMode = 'auto'; fig_pos = fig.PaperPosition; fig.PaperSize = [fig_pos(3) fig_pos(4)]; saveas(gcf,'case1_video', 'pdf') end case 2 % simulate from optimal - needs 40 sec settling time izero = find(deflector.beta==0); % get zero beta [CDmin,imin] = min(deflector.CD(:,izero));% get optimal point deltaHmin = deflector.deltaH(imin); tfinal = 170; % simulation run time (s) thetahat0 = deltaHmin; % theta initial condition - start at optimal VTref = 24; % reference speed m/s VT0 = VTref; % initial speed; fM0 = m*g*(Crr+Crv*VT0^2 ); % initial mechanical force fD0 = VT0^2*rho*CDmin*A/2; % initial drag force Ki0 = (fD0+fM0)/Th2fE; % initial PI control output VTrefchange1 = -5; % reference speed change 1 m/s tVTrefchange1 = 40; % reference speed change time 1 s VTrefchange2 = 8; % reference speed change 2 m/s tVTrefchange2 = 120; % reference speed change time 2 s VW0 = 0; % no wind headingrate = 0; % no heading alphachange = 0; % no gradient talphachange = 0; % gradient change phichange = 0; % no wind heading change tphichange = 0; % no wind heading change % simulate sim('Deflector_esc_UKACC2018'); %% plot results figure(21) plot(deltaH(201:end,1)-20,deltaH(201:end,2),'-',... [0,tfinal-20],[min(deflector.deltaH),min(deflector.deltaH)],'-.',... [0,tfinal-20],[max(deflector.deltaH),max(deflector.deltaH)],'-.',... [0,tfinal-20],[deltaHmin,deltaHmin],... 'linewidth',1.2) h99=gca; h99.YLim=[0.8500 1.2000]; h99.XTick=[0:30:150]; grid on ylabel('$\delta_H$','Interpreter' ,'Latex','fontsize',fs) xlabel('$t$ (s)','Interpreter' ,'Latex','fontsize',fs); % g99=legend('$\theta(t)$','$\delta_H^0$ ' ); % set(g99,'Interpreter' ,'Latex','fontsize',12,'Location','southeast') if(1) % print fig = gcf; fig.PaperPositionMode = 'auto'; fig_pos = fig.PaperPosition; fig.PaperSize = [fig_pos(3) fig_pos(4)]; saveas(gcf,'case2_dH', 'pdf') end % plot results (DETAIL) figure(22) plot(deltaH(201:end,1)-20,deltaH(201:end,2),'-',... [0,tfinal-20],[deltaHmin,deltaHmin],... 'linewidth',1.2) grid on h99=gca; h99.XTick=[0:30:150]; grid on ylabel('$\delta_H$','Interpreter' ,'Latex','fontsize',fs) xlabel('$t$ (s)','Interpreter' ,'Latex','fontsize',fs); % g99=legend('$\theta(t)$','$\delta_H^0$ ' ); % set(g99,'Interpreter' ,'Latex','fontsize',12,'Location','southeast') if(1) % print fig = gcf; fig.PaperPositionMode = 'auto'; fig_pos = fig.PaperPosition; fig.PaperSize = [fig_pos(3) fig_pos(4)]; saveas(gcf,'case2_dH_detail', 'pdf') end % plot results figure(23) plot(VT(201:end,1)-20,VT(201:end,2),'-',... VT_ref(201:end,1)-20,VT_ref(201:end,2),'--',... 'linewidth',1.2) h99=gca; h99.XTick=[0:30:150]; grid on ylabel('$V_T$ (m/s)','Interpreter' ,'Latex','fontsize',fs) xlabel('$t$ (s)','Interpreter' ,'Latex','fontsize',fs); % g99=legend('$\theta(t)$','$\delta_H^0$ ' ); % set(g99,'Interpreter' ,'Latex','fontsize',12,'Location','southeast') if(1) % print fig = gcf; fig.PaperPositionMode = 'auto'; fig_pos = fig.PaperPosition; fig.PaperSize = [fig_pos(3) fig_pos(4)]; saveas(gcf,'case2_VT', 'pdf') end case 3 % simulate from optimal izero = find(deflector.beta==0); % get zero beta [CDmin,imin] = min(deflector.CD(:,izero));% get optimal point deltaHmin = deflector.deltaH(imin); tfinal = 80; % simulation run time (s) thetahat0 = deltaHmin; % theta initial condition - start at optimal VTref = 24; % reference speed m/s VT0 = VTref; % initial speed; VW0 = 20; % wind speed m/s phi0 = 0; % initial wind heading angle VR02 =( VT0^2 + VW0^2 + 2*VT0*VW0*cos(phi0) ); % initial wind incident angle fD0 = VR02*rho*CDmin*A/2; % initial drag force fM0 = m*g*(Crr+Crv*VT0^2 ); % initial mechanical force Ki0 = (fD0+fM0)/Th2fE; % initial PI control output VTrefchange1 = 0; % reference speed change 1 m/s tVTrefchange1 = 0; % reference speed change time 1 s VTrefchange2 = 0; % reference speed change 2 m/s tVTrefchange2 = 0; % reference speed change time 2 s alphachange = 0; % no gradient talphachange = 0; % gradient change phichange = pi/2; % 90 deg wind heading change tphichange = 20; % wind heading change VRchange2 =( VT0^2 + VW0^2 + 2*VT0*VW0*cos(phichange) ); % initial wind incident angle betachange=asin(VW0/sqrt(VRchange2)*sin(phichange)); % wind incident angle %alphachange = 3*pi/180 % 3 deg gradient % simulate sim('Deflector_esc_UKACC2018'); % plot results figure(31) deflector.beta==0 [CDmin(ii),imin] = min(deflector.CD(:,izero));% get optimal point end