Building a Matlab/Simulink Model of an Aircraft: the Research Civil Aircraft Model (RCAM)

hello everyone and welcome to another video today I'd like to continue our discussion on building a nonlinear six degree of freedom aircraft simulator so if you recall in our previous video we discussed the research civil aircraft model also known as the Archaea model we built up the equation / models describing the aerodynamics propulsion gravity etc that comprised the equations of motion for this aircraft this link the theoretical foundation for the aircraft model and now we're actually ready to spend that in mathematical capital and implement in software to generate an honest-to-goodness flight simulator since we spent almost two hours in the previous discussion going over the model in detail we're gonna move at a somewhat accelerated pace today since all we're really doing is translating these equations into lines of code obviously if you haven't seen that video first you might feel more confused than a chameleon and a bag full of skittles so please make sure you watch the videos in order assuming that you've watched the previous video and you're comfortable with the content so why don't we jump over the white board and get started all right so I've rewritten our Flat Earth equations of motion for a rigid six degree of freedom object and here they are right and then last video we discussed what made these equations specific to the Archaea model namely it was how the external forces come in we saw that they were a function of gravity the propulsion / engines and the aerodynamics and also how the external moments were specific to our particular aircraft again the engine / propulsion caused a moment as well as the aerodynamics caused a moment as well as the last couple terms are the mass was aircraft specific as well as the inertia matrix and it's inverse here okay so that was the last video so now what we need to do is we need to go ahead and figure out how to implement all of this in MATLAB remember we spent two hours discussing how to build up all of these functions of what the relationships between the aerodynamics and the state vector were and all that kind of stuff now we're in the position to package this up into a MATLAB function so what's going to do that this the whole purpose of all these complex 9 differential equations is to basically calculate the state derivative dot right so again remember here's our state vector right it said it was uvw PQ r v theta sine well although these Nation ships gives you UV w PQ r dot Phi dot theta dot sine dot right so the left hand side of this what we would want is coming out of this is going to be X dot right okay and if you think about this long enough we now have to figure out what is necessary to compute this basically we want to come up with a function which is basically our our cam model right where its sole job in life is to compute the state dot well what is it going to need to do that well obviously if you look at these equations we see that okay it's always easy getting it you're gonna need constants right we remember that things like the aerodynamic the tables it had a bunch of different slopes and angle of attack where it switches from linear to nonlinear there are all these other constants associated with it so we can bake that in here but then there were things like here this term right here that was just uvw this term right here this is X four five six right so we see that the state vector is also needed so I need to know the state vector at the current time in order to compute the state vector dot okay the other thing that we need it might not seem explicitly written here but if you remember back in our discussion we said things like like the aerodynamic forces right if you were to sit here and adjust the aileron the tailor of the rudder that's gonna affect like the lift drag or the side force and it is also going to affect the the moment here all right and it's all here like these other two if you look at these two control surfaces right the throttles if you crank up these throttles it's definitely going to change this term right here right this is the net force due to propulsion effects and here's a net moment so that's gonna weasel its way in there so long story short you also need the current control right so if I have the state vector at a given time I have the control vector at a given time I have all these constants I should be able to now package all this up into a MATLAB function and have it compute X dot that's our game plan so what I want to do is go ahead and let's go over to the to the to the to MATLAB and basically build this function let's call it actually how about our kam underscore model dot M all right so this thing's like I said its sole job in life it is a function where you pass in the state vector and the control vector and we it will then compute X dot okay and the way we're gonna get that is we're gonna basically go through that exact same ten steps that we derived in the previous video where we built up things like the dimensionless aerodynamic coefficient in F s and then we did all the rotations and yadda yadda yadda in fact I mean just to refresh your memory here maybe I'll flash I'm on the board right here here are the ten steps that we talked about in painful detail last time now what we have to do is translate those equations into actual lines of code so if that sounds good let's jump over to MATLAB and do that alright so we've just got a blank MATLAB script so let's go ahead and define what this function should look like so we said this thing was gonna be called arc am a model and the interface was again we pass in a nine element state vector X a five element control vector u and this thing kicks out a nine element vector which is just X dot so that's pretty good let's go ahead and save this as yep our cam model dot M there we go and now we can just start walking through our ten steps and actually before we do that maybe what we should do is let's have a to small sections to start where we are gonna do a couple of extra things namely and in order to save time I'm not going to type everything in I've actually got this prepared off screen so I'll just copy and paste in the first thing I want to do so the first thing we're gonna do is it might be nice if we break out all the states x1 through 9 into individual variables called x1 through nine just so that we can get away without not having to do all this indexing down the road and again we remind ourselves here that this is the state vector it's UVW p q r v theta sy right in that specific order similarly for you okay that's the first thing that we want to do the next thing that I want to define in our function is all of those constants there were a ton of constants with the Archaea model so again I'm gonna grab from off-screen and I will just go ahead and paste in the appropriate values let's scroll up to the top and maybe just talk through these so um here the first thing the mass of the aircraft the nominal aircraft is 120 kilograms you might be really tempted to define the inertia matrix here and you totally can but I'm gonna defer that as it's gonna be relevant to our discussion a little bit later so tell you what instead of defining that three by three inertia matrix right here you can go ahead and we're gonna do this somewhere else all right let's keep going all these other parameters like C bar this mean aerodynamic chord all of this other stuff everything is pretty standard the only thing that I'll maybe call out is this zcg value um the our cam document has a little bit of an inconsistency in fact I'll flash it up right here depending on if you're looking at table 2.4 or table 2.5 you might see different nominal values of what the our cam document suggests you should use I'm gonna go with the 0.1 C bar value as it seems a little bit more correct I guess and there's less typos associated with this value so anyway other than that there's nothing very fancy about this section right all we're doing here is defining all of these constants like where are the engine's located right these are these application points a thrust of engine one and engine two and then all of our other constants here's our constantly assumed air density remember so we talked about this in the previous video we're gonna assume that the air density is changes with altitude are negligible we're just gonna hold it constant at the sea level value and then again here's all these other numbers the only thing that maybe is also relevant to point out is again remember we talked about in the last video but the a0 number in the our cam document yields a discontinuous lift curl lift the CL alpha slope so we're gonna use this slightly modified value okay so with those two kind of pre sections done we can really just start jumping into the 10 steps that we talked about earlier so the first thing that we needed was uh control saturations so again to save time let me copy this from off-screen and I'm just gonna paste it in here and then we'll walk our way through it so section 1 right control limits and saturation so remember the art camp document said that the ailerons here which was you one can only go between minus 25 and 25 degrees the tail can go negative 25 to positive 10 the rudder can go from minus 30 to 30 and then the engines again this is we're gonna revisit this in a little bit but right now we said that the engine Zeke you actually can't turn them off right they said here's the minimum thrust they provide and here's the maximum thrust excuse me not the actual thrust but this is the minimum deflection the control minimum throttle setting and the maximum throttle setting for engines 4 and 5 ok so then once we've defined them the limits the upper and lower limits the next thing I'm going to do in this next section is really I'm just gonna run each one of the control elements through basically an if check to saturate to see if we're outside of those bounds so for example right if the aileron is deflected beyond or at least the past n value right this is the value that is passed into the function if that is above the maximum value we're actually just going to saturate that value at the upper limit if it's below the minimum value will saturate it at the lower limit so actually after line 81 we guarantee that this value u1 is within the desired control or within the specified control limits for that channel and then we do the same thing for channel 2 3 4 and finally for control number 5 okay so that's section number 1 um what was our next section I think the next section was all of those intermediate variables that we discussed okay so again let me just paste it in so we can go a little bit quicker here so here's the intermediate variables again we talked all about what these were in the previous lecture so I'm not even going to talk too much about them other than again call out maybe some gotcha moments again notice we use 8 and 2 for the angle of attack we do not use a multi quadrant inverse sine so here you may have to think about this a little bit but we're gonna hope that beta is never in ridiculous ranges so this is probably okay um dynamic pressure and then a couple of these turn which are basically just portions of the state vector stacked up okay there we go that's done with number three okay section four was the big one and I think we spent the most time in the lecture last time talking about that but it's actually not a ton of actual code so if we paste this in we can talk about section three so section three um if you notice again all we're doing is implementing all of those equations we talked about earlier so here is our nonlinear lift curve where it is linear until the angle of attack alpha which we have just defined here in line 112 if that is above this this I think it was 14.5 degrees if it's above that we will now use the this this polynomial version and to basically model positive stall right okay that's great and then we just go ahead and do the same thing for our tail right we implement the CL of the tail equations and then we add them together and since the our cam designers were nice and already took into account the difference in normalization normalizing areas we can just directly add them and then again we do the the total drag force here and again remember we had the discussion this is it really dragging the stability in fact the lift drag and cyphers although these are in the stability axis right so ah there we go we calculate all those dimensionless coefficients and then the next thing that we can do is move on to step four and we are going to now just go ahead and dimensionalize those so again I'll copy from my prepared script off screen and here we go we will stack up all of the coefficients right and we will put the appropriate sign and put them in the appropriate either the force in the X direction the Y direction and the Z direction here and we will then scale them all up so what I've done here is just created the dimensional aerodynamic force and the stability axis and then all we need to do is calculate our rotation matrix here right which rotates from the stability access to the body I'm sorry from the stability frame to the body frame and then I can go ahead and apply that rotation matrix to rotate the aerodynamic force from the stability axis to the body axis great okay so moving on let's move on to step number five step number five was I believe getting the dimensionless aerodynamic moment coefficients about the aerodynamic center in FB oh that's a lot to say right okay so here we go step 5 so again we are just implementing the equations that we discuss in our last video so remember that ADA vector was sort of our static stability that we discussed so I've just broken it up into a couple of intermediate variables but again you guys are all good at MATLAB at this point this is just making that happen I'm gonna define that three by three matrix which we kind of use that abuse of notation and called it DCM DX but it really it's the it's how the moments change with respect to a state namely X four four five and six right that's what this matrix measured and then finally here was our quote-unquote control authority or control effectiveness sort of matrix again this one was going to be multiplied with the controls as you can see right here controls one two and three so again this is just implementing those equations that we talked about in detail last lecture so I'm just going to gloss over them and just show the code that we need to get them entered in okay that was the the tricky part in terms of defining the code because now we've got the dimensionless aerodynamic moment about the aerodynamic center in the body frame so step six is really just scaling dimensionalizing that dimensionless coefficient and again maybe we will call out that the thing that's nice about this is the our camp designers again decided to use the exact same normalizing length for all of the moments so really um yeah you don't really have to do anything here you can just go and normalize all of them by the exact same normalizing area and the exact same normalizing length to get the actual moment here so this is the moment in terms of Newton meters about the aerodynamic Center expressed in FB right and espress in the body frame all right so moving on to step 7 step 7 is just performing the moment transfer to move these moments from the aerodynamic center to the center of gravity and again maybe we'll just remand call out that I believe the our camp document has a sign error in their expression so we are going to use our derived expression to do this moment transfer and maybe I'll also call out that MATLAB has a nice cross function that you can use to do the cross-product so uh nothing much to say about this since we are again we already spoke quite a bit about this in our previous video so let's move on to step number eight which again I'm prepared off-screen I'll just paste in the code and we'll talk about it so here we go here are the engine forces and again maybe the thing to call out is that the force produced by engine one and two right we see that they are directly related to the through the input throttle setting u4 and u5 now remember way up at the very top we had saturated this to a minimum and a maximum value and the thing that was interesting about this as the minimum value is above zero so we see that in this current implementation of the our cam model there's no way to stop the forces generated by these two engines right they're always going to be providing some minimal amount of thrust um and again I'll call that out because we'll come back to that in a little bit but otherwise assuming that is all set we can go ahead and calculate the actual force that those two engines are producing as a vector add them together to get the total propulsive force Express in the body frame and then use our cross-product again along with the difference between the position vector of the application point of thrust and the CG of the aircraft and we'll get the moments that those two engines produce about the body about the center of gravity Express in the body frame ok so there we go that was step eight now step nine which was our gravity effects was really easy again we talked about this already so there's no need to discuss here uh what else let's keep going our last step now this is where it gets sort of interesting I'll go ahead and copy it from off-screen and paste it on here to take a look so let's go ahead and chat about this so the first thing you note and maybe we'll call this out now is that right here is where I define the inertia matrix again this is a matter of personal preference you could have defined this inertia matrix way up at the top where we were defining all the other constants like mass what not the only reason I did it down here because I think there's a little trick that's worth mentioning and I think it helps visualize why we need this trick so if you remember down here on line 239 this is basically implementing Newton's second law of rotation for our Flat Earth equations of motion right and notice that we need the inertia matrix and we also need its inverse of the inertia matrix if you notice now what I've done here on line two to nine here is I've actually gone ahead and off line I have Hart I've computed the inverse of this constant matrix and hard-coded it in here right um you could alternatively say inv I begets inverse of IB right this is another way to do this right this is another way to implement line 229 where we're just getting the inverse now the only reason I'm calling this out right now is this is kind of a stupid thing to do right because the inertia matrix is constant right it's never ever changing so if we were gonna call this our cam model dot M function like thousands of times say within a simulation why would you waste time asking MATLAB to compute a 3x3 inverse every single time and it's just gonna give you the exact same answer every single time so instead to save that computational time instead of asking my lab to compute the inverse I'm just gonna compute it offline one time and then hard-code it in here to save me some effort so the other thing that maybe note is I don't think we discussed the structure of the inertia matrix in our previous video and now might be a good time to do so if you notice here this matrix is symmetrical and there are a lot of off diagonal zero terms this is consistent for an inertia matrix of a vehicle or any object which is symmetric about the XZ body plane which it totally is for an aircraft or at least for most aircraft right if you split an aircraft right down the center the starboard side looks exactly like the port side it's just a mirror image reflection the other thing that maybe we'll call out right now is that the Arkham document again seems to have a bit of an inconsistency if you look at how they define the inertia matrix on page 12 versus how they define the inertia matrix on page 91 you'll see that there's an inconsistency in the thousands and ten thousands place of one of these entries and furthermore there's a further type on the the entry that has an inconsistency it is the X of the I XZ entry and you see there's a cite typo there so due to those two typos on page 91 I'm going with the version that's on page 12 here which is what we've implemented okay so enough talking about the inertia matrix again the way the inertia matrix is used is right down here on line 239 where we basically implement Newton's second law to calculate X 4 5 & 6 . similarly on line 235 we have Newton's second law of translation which gets us X 1 2 & 3 dots so um lastly down here on line 246 we are implementing Euler's kinematical equations which again we had in previous video discussing orders kinematical equations and this is just our choice here to make our lives easier and make the code easier we could have had a whole other section on implementing Pisans kinematical equations or the quaternion kinematical equations but to make our life easy for this first cut let's go ahead and use Euler's kinematical equations so at this point we have again this is X 1 2 3 . this is a 3 by 1 vector this is X 4 5 6 dot which is again a 3 by 1 vector and then finally this is X 7 8 9 dot which is another 3 by 1 vector so all we need to do is stack them one on top of each other to get this 9 by 1 vector which is the state derivative dot so that's all there is to it if we come now up to the very top we see that you pass in the current state vector in the current control vector and in just 250 lines of code we can compute the state derivative dot of this full nonlinear six degree of freedom twin-engine 120 kilogram commercial transport jet so let's go ahead and save this and we've got the our cam underscore model dot M function which gives us this functionality so what do we want to do with this now all right so we've fully defined this MATLAB function right our cam model dot but this is not actually the entire aircraft right because if you look at this right now X dot and X these two signals are not related the way that we've drawn this right and we know that these two have to be related right they're related through one integration so in other words if I go ahead and now stick an integrator right here what this is going to now output is X right you integrate X dot you get X so in order to make this consistent I need to actually wrap this back like such right so now this overall system has the control surfaces coming in it's going to calculate the state derivative dots it's then going to integrate it and feed that back and use that for the same for that for the calculation at the next time step right so this entire picture right this block diagram is consistent and this is our aircraft so if I were to just outline this maybe in blue like such and maybe what I'll do is let me pull this vector a little bit further out all right so now this thing in blue is the non-linear six degree of freedom aircraft simulation right because what comes in here is on the left side you have the pilot inputting the whatever the aileron the tail the rudder the throttles whatever those are and then all this integration is going to occur and it's going to basically compute the state right of the vehicle how fast it's going its angular velocities and its oiler angles and that will come out on the far side so this right here is our model of our aircraft right this is our our cam aircraft so what we need to do is let's go ahead and see if we can integrate this into MATLAB Simulink and build a block which implements this in Simulink all right so let's go ahead and fire up simulate and start ourselves a blank model and the first thing we want to do with this model is let's go ahead and save this in the same directory as the Archaea model dot M file so maybe let's call this our cams simulation hit save okay and we can now start building the diagram that we had on the board so I'm gonna start with the integrator that's what I like to do and what's coming into the integrator is our signal X dot which maybe let's label this explicitly this would be our nine by one vector we integrate that and coming out of this should be X which is also a nine by one vector the last thing that we should maybe talk about before we leave the integrator is the initial condition here let's just set this to a variable called X zero and I will come over to MATLAB and let's start ourselves a script whose sole job in life it is is to initialize constants whoops initialize constants for the our cam simulation right so let's start off with our good old friend clear CLC close all and then let's go ahead and define constants here in this first cell so x0 is going to be again I've got this prepared off screen so you don't have to watch me typing things in this is the initial condition so remember this is UV w P Q R Phi Theta sighs so we're going with a u velocity of about eighty five meters a second which is about 165 knots and we're pitched up at about five point seven degrees so this is the initial condition okay coming back to the Simulink model we've got that to find the other thing that we need is probably our control vector so we need to know what control surface deflections are we trying to fly at so I'm gonna come over on the far left and let's just set them to be a constant for now so this constant let's set this just again to another variable of let's just call it U and we'll come back to our script where we are initializing constants and we'll pick our initial you so the initial you that I'd like to use is this so remember this is basically saying zero ailerons the tail is actually down at about negative five point seven degrees there's zero rudder and then both controls for the throttle are they're set symmetrically at point zero eight and if you recall earlier the limitations on the throttle the minimum value was 0.5 times PI over 180 which evaluates to about 0.009 so we're not saturating the controls in any of the channels and both of these throttles again are applied symmetrically so both engines have the exact same throttle setting okay great so we've got that set up so coming out of this block is going to be the signal U which again is a five by one vector so at this point we basically have everything we need if we branch off this signal X and bring it back you'll see that now what I need right here is I need to be able to call the our cam model dot M to count to calculate X dot right so that's what we need there's a couple of ways you can do this one of them is you come here to your Simulink library browser come here to user-defined functions and grab this MATLAB function and drag drag and drop them into the model and if you double click on that you'll see what this MATLAB function block allows you to do is basically write MATLAB code and input that so if you come over here to the our km model dot M file what you can do is control a select everything copy it and come back to the MATLAB function block and paste all of that code inside and now if you save this and you come back to the Simulink model you'll see what that does is it basically allows us to perform that calculation this is the same calculation that our cam underscore model dot M was making so one way to skin this cat is you can just drop this in right like this so let me hook everything up just for completeness and then X comes over here whoops let's bring that up okay like that oops sorry there like that let's hit save um so this is one way you can do this this should work right now the only problem though the thing that I don't like about this is if you think about this long enough unfortunately right now we have two versions of the code running around so there's this version right here so this one on the right this is the code that is inside the MATLAB function block and it's completely separate from the code that's in our cam underscore model dot M both of them have the same contents but they're different instances well for those of you who are software developers I think this is probably throwing up a bunch of red flags in your head because I don't know about you but I detest having the same copied in two separate locations right it makes maintainability and nightmare like for example what if we came here in the our kam underscore model dot m and somebody said you know what the vehicle doesn't weigh 120 thousand kilograms that's not the mass it's actually 110 so you come here in our kam underscore model dot m and you update the mass to be the appropriate value and you hit save and a whoopsie-daisy look at this it's not automatically updated in the Simulink diagram so basically this Simulink model is now inaccurate it's not reflecting the current value that you want so again this is not a great idea to have two instances of the same code floating around so instead i really don't want to use this block so let's go and close this switch the the mass back to the nominal value and come back to the simulate mall and delete this so really i don't want to use this instead what i want to do is i really want to have matlab call this particular our camp underscore model dot m function and use it to calculate the state derivatives rather than copying it into a matlab function block so the way we're gonna do that is let's go again come to here to user-defined functions and grab an interpreted matlab function and drag and drop now if you look at the block parameters for an interpreted matlab function you see what it allows you to do is it allows you to call any matlab function like sine cosine tangent whatever or check this out you can also call user-defined functions right as long as this on the mountain a path and matlab knows where this function is it can call it so what we can say right here is we could say all right I want you to call our kam underscore model and then here's the one issue I need to pass in parameters so the thing that makes is interesting is if you look at this the interpreted matlab function only has one input nubbin coming in i need to pass information about both the state and the control vector to our kam model dot m alright so the way to get around this and the way to do this this is an interesting little trick is let's go ahead and grab ourselves a MUX or a multiplexer okay there's multiplexer what I'm gonna do is I'm going to now use it to stack the state vector on top of the control vector like this okay so maybe let's double click this to get the label over here so we can see what's going on so you've got X on top and you on the bottom so coming out of this thing is actually this giant mega vector called X U which is now a 14 by 1 vector right so what we can do now is if i hook this up like this okay this is where we're almost there I need to open this up and take a gander at the input parameters now here's where we have a slight naming collision if you're reading the help documentation on this interpreted MATLAB function notice here that within the scope of the interpretive MATLAB function you right here if I type in the the symbol u lowercase u this refers to the inputs to the block or the input signal coming into this interpreted MATLAB function block which we see in our case it's actually a 14 by 1 element vector now the problem here is our function our camp underscore model if you look at it up here in the MATLAB script right it expects to separate inputs it wants a 9 element vector first followed by a 5 element vector so you we have two input two separate vectors so the way to get that to happen here in the block parameters and the interpreted MATLAB function is you say u 1 2 9 and then comment u 10 to 14 so hopefully if you stare at this long enough it's gonna make sense right what this does is it says ok take the first elements of the 14 element input which in our case is the states and pass that as the first argument to our cam model then take inputs 10 through 14 which you see right here is actually the control vector U and pass that into our cam model as the second input so again I think I'm I think I'm rambling at this point I think hopefully this makes sense if you if you look at how these are lined up right okay so with that being said I'm gonna hit apply ok and I think what's I think we're ready to go so let me save the art camp simulation let's come back here do our initialization script I guess we should have saved this I guess let's save this here in the same directory and call this thing like initialize constants okay and oh maybe let's do one more constant let's say T a TF the final simulation time let's make this like 60 seconds and then I'll come back to the Simulink model and change this instead of ten seconds let's make this TF okay so we simulate for 60 seconds okay um what else should we probably do let's save some of the data so for example let's save the state and the control vector so what I'm gonna do is let's come here to our sinks and let's get a two workspace block and drop this in here and I'm gonna call this thing sim X cuz that's the simulated version of the state vector and then over here I'll hook this up call this thing sim you okay and it will use a time series and I'll go ahead and hit save and now if I come over to my initialization script I can go ahead and now run the model from the script using the sim command so what we need to do is we need to say sim um what did we call that thing was the name of our simulation oh here it was called our cam simulation dot Philips simulation dot SLX there we go and if I now run this script it should define the constants run the model and generate the data and there we go we've got SIMEX and sim you created maybe now is a good time to maybe mention I'm running MATLAB 2018 B I think if you're on a newer version of like MATLAB 2019 or 2020 or anything newer these variables may be saved in a different data structure I believe like an out data structure so if you're on a newer version of Simulink you may need to modify the code slightly but I'll assume you're you're able to do that so now that we've got the data what we can do is let's go ahead and plot the results so let's go ahead and again to save the time I'm just gonna go ahead and copy and paste from off-screen so again all this is doing is just unpacking all of the data from these two time series that were created and then let's go ahead and plot all of this again to save ourselves time I've got this all prepared so for example I'll plot the first control input this should be the aileron and then let's go ahead and do this same thing for all the other control inputs on the same figure so here we go so here's plotting aileron tail rudder control one control two and then let's also plot how about these states because that's what's more interesting so if we go ahead and do that again I've got all this prepared off screens so that you're not gonna have to just sit me watch me typing in a whole bunch of stuff I'll paste this in and we'll whoops learned it ago okay so we'll take a look at this so again um yeah so in the subplot I've got the first three as being X you know X 1 2 3 which are UVW and then PQ R and then v theta side ok so if we run this this should now do everything it should initialize the model run the model and then plot up the results so here we go let's take a look uh whoops what happened here oops that should be a 5 sorry hold on let me let me update that and that should be u5 sorry ok let me run this again I had a small little bug there ok it's running and let's take a look so here's the first one here's all the controls as expected since we just had a bunch of constants right you remember look here this is just a bunch of constants you can see that ok no aileron we had the tail at the right spot we had no rudder and we had a constant throttle on both now let's see what happens what does the vehicle do now let's pull this over and take a look so here's X 1 X 2 X 3 so here's U V W so what's interesting is we started at 85 meters a second but actually we lost some speed then we increase we gain some speed all the way up to G's a hundred and thirty five ish you know what is this you know uh come on flash up there ah keeps popping up and down when I hover over ya 136 m/s and it looks like it was stabilizing there's no motion in the lateral directional states which kind of makes sense right we started out with the aircraft only having longitudinal motion so there's no W component and there's no P and there's no R and there's no phi and no sigh which again makes sense so again it looks like the vehicle started slow and it did a little bit of oscillations and it looks like it's about to probably come to hopefully to steady-state so for example if we ran this thing longer so for example let's come back to our MATLAB script let's run this instead of 60 seconds maybe let's triple the time let's go for 3 minutes of simulation and run this again and hopefully we'll see what happens now and again I'll pull the state trace over and yeah this looks good so it looks like at that throttle setting it looks like the vehicle will eventually come to steady state here you know I'm around 120 meters a second and yeah this is very interesting right now what makes this kind of cool is we can now go ahead and start modifying this too to simulate different environments so maybe the first thing we should do is to make this perfectly clear let's make a subsystem out of this so remember all these blocks this was the aircraft so I'm gonna select them all right-click and then let's go ahead and create a subsystem so this is basically the nonlinear six degree of freedom our cam aircraft plant model right so you stick in some control surface deflections and out comes the state so now we can start playing some games so for example let's add like maybe an aileron pulse and see what happens so what I mean by that is tell you what let's let's just do this quickly this should be pretty easy and kind of fun let's go ahead and break this line here let's make a summing block and we will have the you coming in here and then let's maybe make some small Delta U I'm gonna make this whoops I want actually like a step let's grab a step real fast that should be under sources and let's grab a step and I don't know let's let's make a step at six at 30 seconds how about so step a 30 seconds let's step - maybe like a 5 degree aileron and then let's hold that for just you know two seconds so maybe 30 plus so 32 seconds and then we'll go to negative 5 degrees and if we now add those together okay that should be a pulse in aileron and now we don't want to pulse anything else so let me go ahead and make a get a MUX here we'll have five inputs and what we'll do is here's the pulse in aileron and then the pulse in all the other channels let's just make them a whole bunch of zeros so grab a constant put a zero and we're going to basically do something like this just put zeros into all these whoops come on zero zero zero zero okay and connect these up okay so this signal coming out of this thing let's maybe turn this into a subsystem this is like a delta u alright so if i hook this in what we should be doing now is we're gonna now fly along you know in the longitudinal motion for just 30 seconds and then at 30 seconds we are gonna crank the ailerons to 5 degrees hold them for two seconds and then bring them back let's see what this does to the simulation now so if I again hook all this up hit save and now we can come back to our script and rerun this puppy and give it a second and let's go ahead and pull over the controls first so here's our plot of all the controls so as you can see we're flying along for 30 seconds and then we pulse the ailerons you know a 5 degree aileron deflection for two seconds and then we bring them back to zero okay now what's interesting about this is if you look at the state trace information again there's no longitude lateral directional motion for 30 seconds but then after you kick the aileron look at this you start inducing some sideways velocity you start inducing a roll I mean rate that's exactly what the ailerons are supposed to do right they induce a roll rate therefore it also induces a yaw rate kind of as a second or secondary after effect and look at this this is interesting so luckily for us it looks like even with this aileron kick um the system becomes stable and stabilizes out however look at this it stabilizes and now it's no longer heading north any longer right psy is moved to this you know almost negative 0.3 radians so this is really interesting there's a lot of interesting behavior now in this model that we're gonna explore later but long story short at this point we've set up a pretty good representation of an aircraft here in this block so um maybe what we should do right now is revisit that discussion of these limits and the saturation all right so here's the issue with the way the limitations and the saturation on the control vector is currently being implemented let's say you wanted to do something like simulate an engine out or something like that so an engine out would be let's come back to our defined constants we could maybe say let's kill engine one so zero control input so I want to basically pretend like engine one provides no thrust and no force on no moment on the vehicle so I could do this however the problem right now is if we trace that signal in the our cam model in our Simulink model right so coming out of this u4 is zero u4 still zero u 4 0 coming into everything u 4 is zero u 4 is still zero now coming into the interpreted matlab function we're gonna pass into here right u4 as zero so now we got to come over to our cam model so here's our our cam model right down here u 4 is still 0 so this thing right here on line 18 that's gonna be 0 right now the problem is it's 0 here which is which is good for our engine now but however when you come down here to the section where we introduced control limits and saturation right here we're gonna trigger this statement right so we're gonna land here on line 98 and u4 is not going to be set to 0 u 4 will go to minimum so basically in this implementation there's no way for us to shut off the the force or the throttle to one of these engines completely and I think it's a little bit insidious in the sense that this saturation limit is buried all the way inside the our cam underscore model dot M I mean as long as you're aware of it maybe that's fine so what I suggest instead is in instead of having these limitations or these trol vector saturations buried inside the Arkham model dot M function maybe what we can do here is actually let's just copy all of this code okay which is our saturations I will just copy this and maybe let's comment this out inside the model here so I will comment this out comment okay so now our cam underscore model dot M is not doing any of the saturation instead what we could do is let's paste some of this code into our initialized constants so in here where we're defining all the constants maybe I'll put in another little small subsection where we're gonna say um define minimum / maximum control values so I'll paste some of that in here what I want to keep is the top stuff I don't want to do any of this all this if statements right this was how we were trying to enforce those upper and lower limits so I'm gonna delete all that all I care about is I want to keep these definitions on what is the maximum and one of the minimum values of all of these channels view one two three four five okay so if I do this and now come over to the Simulink model what you can do is inside here on the you signal what I can do is let's go ahead and first break this out so I'm gonna grab a D mocs v element D MUX okay let's do it like such and maybe I will disconnect this temporarily hook this up like so and now what I'm gonna do is let me let me let me blow this up so we have a little more real estate to work with I'll move this up okay okay what I'm gonna do is I would like to now implement saturation here so I'm gonna say saturation I'm gonna grab a saturation block saturation does exactly what I want to here this is going to be u1 and what I want coming out of this guy is basically u1 saturated which is respecting and obeying those limits so now if I come into here I can say this upper and lower limit what did we actually call this let's come back to the script here u 1 min and u 1 max so what I want can do is in this saturation block the lower limit is u 1 min the upper limit is u 1 max right so that will that will implement saturation let's do the same thing for all these other channels right and we'll hook them up and then all we have to do is change these values in here to the appropriate one so now this is u 2 & u 2 and then this guy is u 3 whoops come on ah geez u 3 & u 3 and then finally here's the engines right this is now you four and u 4 and then finally we're at their last one u 5 + u 5 ok and now let's MUX these all up together and again this is now five elements and if we hook this up hopefully what we've got now is if we so again select all these let's get them all connected oh come on don't tell me it's not gonna be smart enough about that all right fine we're gonna have to manually hook all these up but I think you get the picture at this point and I guess these are not u1 sets anymore right this is u2 saturated u3 saturated you for saturated and u5 saturated so now if we do something like this and hook this up right we've basically moved the saturation and the limiting out of the are cam underscore modeled on em and we've done it explicitly here where it's very easy to see what's going on and more importantly what we can do now is in our initialization script we can change these to appropriate value so example if we wanted to be able to shut off you for completely we need to be able to change this minimum here to instead of being 0.5 this needs to be able to go all drop all the way down to zero alright so now we're allowing the engine for to be completely shut off so that when we go ahead and put in this value it will turn itself off so now if I run this simulation again well you know I don't anticipate this going well because we're flying for three minutes with one engine well it ran let's go ahead and take a look at the results so here you go here's our aileron pulse and then look at this here u4 is sitting here at exactly zero right so the engine for is complete lay off and hopefully what we've done properly is since we we handle the saturations all properly the force from the engine that's actually acting on the vehicle is actually zero and as you can see we're in big trouble here when that happens look this is like yeah stuff is going kind of wrong although surprisingly this is not too horribly wrong it does almost look like we actually end up in a steady state Wow so I'm actually shocked I guess at this low throttle setting it looks like we are actually flying stable II although if you look at sigh it looks like we're just in a turn right we are in a counterclockwise westerly turning you know we're basically going around in circles at this point but again that makes sense we've got one engine on and one engine off so this is um this is great I think I think we're at a really good spot maybe what we should do is I will whet your appetite a little bit by you know looking at this graph you know we're trying to interpret what is going on right what is the aircraft doing it's a little hard at least in my opinion to just stare at state traces and try to understand what the vehicle is doing so what we're gonna look at in a future video is how can we take this Simulink model that we built with our aircraft and maybe instead of just sending data out to be plotted can we send that out to some kind of a visualization package to visualize like a flight simulator to visualize what the vehicle is doing and how its responding and even further down the road if instead of just having inputs right here of just constants or doublets or pulses or sine waves or things like that wouldn't it be awesome if we could just put in like read from a joystick or something like that so we could actually have a pilot hopefully fly the vehicle by in pointing joystick commands our system the simulation the flight simulator that we developed will be doing all of the engineering calculations in computing the states and then we will send that information out to a flight simulator to visualize it so really if you think about this we are gonna be building our own honest-to-goodness full non-linear six degree of freedom our camp simulator slash flight simulator right so this is gonna be pretty fun coming on down the pike so stick around I'm very excited for that so with that being said I think this is a good spot to leave it I hope you enjoyed the discussion today and if so I also hope you'll consider subscribing to the channel surprisingly if you just scroll down a little ways and click on that subscribe button it really does help me continue making these videos please also let me know what you think about this video or any other videos down in the comments I really would love to hear from you with that being said I hope to catch you at one of these future discussions and until then I think I'll sign off so talk to you later bye

Channel: Christopher Lum

Views: 14,865

Keywords: Simulink aircraft model, Matlab aircraft model, flight simulator, aircraft model, RCAM, dynamic aircraft model, simulating flight

Id: m5sEln5bWuM

Length: 51min 29sec (3089 seconds)

Published: Wed Apr 29 2020