diff --git a/sw/in_progress/airspeed/ultrasimplesim.m b/sw/in_progress/airspeed/ultrasimplesim.m new file mode 100644 index 0000000000..30e77504d6 --- /dev/null +++ b/sw/in_progress/airspeed/ultrasimplesim.m @@ -0,0 +1,244 @@ +clear +clc +close all + + + +dt = 1/60; +simsteps = 16000; + +timespan=0:1/60:100; + +t = 0; +airspeed = 14; +acceleration = 0; +vclimb = 0; +altitude = 100; +pitch = 0; +alpha = 0; +pitch_sp = 0.2; +throttle = 0.5; +power = 0.25; +elevator = 0; + + +% Airspeed command +airspeed_sp = ones(1,simsteps).* 11; +airspeed_sp(1,6000:12000) = 16; +airspeed_sp(1,13000:14000) = 7; + +% Altitude command +altitude_sp = ones(1,simsteps).* 100; +altitude_sp(1,2000:4000) = 125; +altitude_sp(1,8000:10000) = 125; + +battery_good = ones(1,simsteps); +battery_good(1,15000:16000) = 0; + +% Save state +X = zeros(simsteps,13); + +% Very sofisticated Aircraft Model +Vmax = 26; % Level flight full power +Vmin = 8; % Stall +Mass = 1.0; + +% Controller States +pitch_cumsum = 0; +airspeed_cumsum = 0; + +controller = 1; + +for i=1:simsteps + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Simplified A/C long dynamics with stall and energy + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Pitch loop without integrators + elevator = (pitch_sp - pitch) * 15; + if (elevator > 0.4) + elevator = 0.4; + elseif (elevator < -0.2) + elevator = -0.2; + end + % Follow elevator but also align with the flow + pitch = pitch + elevator * dt - alpha * dt; + % Angle of attack from 0 deg at vmax to 17 deg stall at vmin + if airspeed > Vmax + alpha = 0; + LoverD = 6; + elseif airspeed < Vmin + % from Vmin downto 3/4Vmin ... increase alpha from 17 to 45 deg + lin = (Vmin-airspeed) / (Vmin/4); + if (lin > 2) + lin = 2; + end + alpha = (17 + (45-17) * lin) / 57; % complete stall + LoverD = 8 + (4.1-8) * lin; + else + alpha = (Vmax - airspeed) / (Vmax-Vmin) * 17.0 / 57; + LoverD = 8; + end + % Path angle + gamma = pitch - alpha; + vclimb = sin(gamma) * airspeed; + altitude = altitude + vclimb * dt; + % Using the pitch angle... + P_climb = vclimb * Mass; + % Lift over Drag + P_drag = airspeed / LoverD * Mass; + % Motor Energy + if battery_good(i) == 0 + throttle = 0; + end + power_sp = throttle ^ 2; + power = power + (power_sp - power) * 0.1; + % Full throttle = same energy as vmax glide + P_motor = (Vmax / LoverD) * sqrt(power) * Mass; + % Total energy + P_tot = P_motor - P_drag - P_climb; + + acceleration = P_tot / Mass; + airspeed = airspeed + (acceleration * dt); + + if (airspeed < Vmin/2) + airspeed = Vmin/2; + end + if (altitude < 0) + altitude = 0; + end + + %%%%%%%%%%%%%%%%%%%%%%%%% + % Run Control 1 + + if controller == 1 + Throttle_min = 0.2; + Throttle_max = 0.9; + + Pitch_min = -20/57; + Pitch_max = 30/57; + + + + % Traditional with integrators: + + vclimb_sp = (altitude_sp(i) - altitude); + if (vclimb_sp > 3) + vclimb_sp = 3; + elseif (vclimb_sp < -3) + vclimb_sp = -3; + end + throttle = vclimb_sp * 0.33; + if (throttle > Throttle_max) + throttle = Throttle_max; + elseif (throttle < Throttle_min) + throttle = Throttle_min; + end + % P-gain + pitch_sp = 0.1 + vclimb_sp * 0.05; + if (pitch_sp > Pitch_max) + pitch_sp = Pitch_max; + elseif (pitch_sp < Pitch_min) + pitch_sp = Pitch_min; + end + + %%%%%%%%%%%%%%%%%%%%%%%%% + % Run Control 2 + + elseif controller == 2 + Throttle_min = 0.15; + Throttle_max = 1; + + Pitch_min = 6/57; + Pitch_max = 18/57; + + + + % Traditional with integrators: + + vclimb_sp = (altitude_sp(i) - altitude); + if (vclimb_sp > 3) + vclimb_sp = 3; + elseif (vclimb_sp < -3) + vclimb_sp = -3; + end + % I-gain + pitch_cumsum = pitch_cumsum + vclimb_sp * 0.001; + if (pitch_cumsum > Pitch_max) + pitch_cumsum = Pitch_max; + elseif (pitch_cumsum < Pitch_min) + pitch_cumsum = Pitch_min; + end + % P-gain + pitch_sp = pitch_cumsum + vclimb_sp * 0.05; + + err_airspeed = airspeed_sp(i) - airspeed; + airspeed_cumsum = airspeed_cumsum + err_airspeed * 0.0002; + if (airspeed_cumsum > Throttle_max) + airspeed_cumsum = Throttle_max; + elseif (airspeed_cumsum < Throttle_min) + airspeed_cumsum = Throttle_min; + end + throttle = airspeed_cumsum + err_airspeed * 0.08; + if (throttle > Throttle_max) + throttle = Throttle_max; + elseif (throttle < Throttle_min) + throttle = Throttle_min; + end + + end + + % Save For plotting + X(i,:) = [ airspeed_sp(i) airspeed altitude_sp(i) altitude throttle power pitch_sp pitch alpha gamma elevator vclimb_sp vclimb]; + +end + + + +figure +subplot(6,1,1) +hold on +plot(airspeed_sp,'r'); +plot(X(:,2),'b'); +plot(ones(simsteps,1)*Vmin,'k'); +plot(ones(simsteps,1)*Vmax,'k'); +title('Airspeed') +grid + +subplot(6,1,2) +hold on +plot(X(:,7)*57,'r'); +plot(X(:,8)*57,'b'); +title('Pitch') +grid + +subplot(6,1,3) +hold on +plot(X(:,5),'r'); +title('Throttle') +grid + + +subplot(6,1,4) +hold on +plot(X(:,9)*57,'g'); +plot(ones(simsteps,1)*17,'r'); +plot(X(:,10)*57,'k'); +title('Alpha Gamma (degrees)') +%legend('1','2') +grid + +subplot(6,1,5) +hold on +plot(X(:,12),'r'); +plot(X(:,13),'b'); +title('VClimb') +grid + +subplot(6,1,6) +hold on +plot(altitude_sp,'r'); +plot(X(:,4),'b'); +title('Altitude') +grid +