Numerically Linearizing a Dynamic System

hello everyone and welcome to another video today I'd like to talk about numerically linearizing a dynamic system and maybe let's break this title down a little bit before we really jump into it so hopefully the dynamic system portion of this is fairly familiar to those you who've been watching the video series so far right I have some kind of dynamic system which has an input vector of U so these are all the inputs that go into the system and then the system through some set of dynamics then those u inputs influence or introduce some kind of a change in the states X right now typically this relationship here is nonlinear in the sense that the relationship between the states of the control inputs the most famous and most common form of that the equations that would govern this are basically a set of differential equations right which looks something like this right in a very generic form okay this is also maybe well let's put another one more descriptor on this this is an explicit form of first-order differential equations it's explicit in the sense that you're able to solve for the term X dot okay so we have basically have X dot is equal to some nonlinear function of the states and the controls right so here's my dynamic system and now it's nonlinear because this F term right this could be really anything you know that you're you can imagine write this F could be describing an aircraft a simple pendulum a cart a set of virus dynamics you know could really be anything right it's very generic now what I want to do here let's move back in the title I want to linearize this so I want to be able to through some process called linearization I want to be able to turn this nonlinear kind of complicated model into a somewhat simpler model okay and in fact I want this to be a linear now dynamic system okay so now this thing has inputs of Delta U right which is the perturbations from some trim point and this then through a similar set of linear dynamics would influence some set of changes Delta X right which again represent perturbations from a trim point the relationship here now I want this to be a linear relationship in the sense that I would like to get a set of first order differential equations which look like Delta X dot is equal to a delta X plus B Delta u right so this here is my linear you're familiar linear time-invariant system right an LTI system okay so that's the process right that's what linearization does it says you take some nonlinear dynamic system you find some kind of a trim point and then we're gonna linearize that about that system about that trim point to yield a linear dynamic system right now the last part of this title is this adverbial in the front numerically right so what I want to look at today is this linearization process we've talked about this before normally the way this is done is this is an analytical process in the sense that you basically compute a whole bunch of Jacobian matrices you take the partials of F with respect to X you take the partial of F with respect to you and you go through this rigmarole of linearization by basically analytically calculating a whole bunch of partial derivatives well that's fine and dandy if this F is easy to deal with and if you have an analytical representation I want to look at the case where what if this F is very very complicated or what if this F doesn't even have an analytical representation what if the F is like a black box right what if the F is a bunch of look-up tables cobbled together maybe with some analytical derivative or analytical descriptors as well but if this thing is just some box that you don't know what's going on under the hood can you still take this crazy nonlinear complicated system this black box and turn it into a linear model so this linearization process we really want to put the adjective numerical linearization right so that's the game plan for today that's what this title means of demands take this nonlinear system and through some numerical set of tools how can we get a linear dynamic system okay so we're gonna look at three ways to do this today so method one here we are going to be using a set of numerical techniques to numerically calculate partials and we're gonna see later on down the road I think what that means for those of you who are familiar at control theory if you already have an inkling of what this means but that's the first technique I want to look at okay and then in method two what we would like to look at is actually MATLAB and Simulink have a set of tools to allow this process to occur once we do method one we're going to have a good physical mathematical understanding of what of how you might be able to turn this into an algorithm well now that I've been Simulink i've already done this so we're going to look at MATLAB Simulink and mainly it has a system called the linear analysis tool okay and then method three is there's more than one way to skin this cat and again matlab and simulink have another function called Lin mod which will allow you to also perform this so this is our road map for today I want to spend most of our time looking at method number one where we're going to try to lay the mathematical framework and the foundation for how you might be able to develop an algorithm to do this in fact we will develop an algorithm we'll develop our own function to do this linearization and then we want to compare our results with these two different other methods that are provided natively in MATLAB and Simulink okay so with that being said there's going to be a fair bit of background material that where we've already covered some of these tools so when we get to that portion of the lecture I'll flash up the relevant prerequisite lecture that I hope you've already watched but that being said I just wanted to warn you ahead of time that we're gonna move a little bit quickly in the discussion today because I'm gonna assume that you've sort of been following the lecture series up to this point and you're familiar with some of these but again if you're not what I would recommend is I will pause the video I'll flash up a reference for you to watch and when we get to that point of lecture you can just pause the lecture watch the prerequisite video and then continue on so alright so that sounds like good fun why don't we go ahead and just jump right into it then all right before we jump into this I want to talk a little bit about the implicit form of Don linear Odie's okay and what I mean by that is if you look at our formulation of our nonlinear set of equations of motion here right our nonlinear set of ordinary differential equations right this is something you've probably seen before right where X thought is some non linear function of the states in control right but I called this earlier an explicit form because you are able to explicitly solve for X thoughts so what you are doing you're writing down the equations of motion which describe this system it just so happened that those equations of motion came out in such a way that you could explicitly and nicely solve an isolate X dot and move it over to one side of the equation right but that's not always the case right so in this form maybe we should write is let's just quickly write this down so again you have the normal kind of traditional form that you would usually see and in this case X dot right this is just some vector in RN right there are n states so there are n obviously X dots and X here is obviously also an element vector and U is a M element vector right so there are n states November 4 n 4 November States M mike-mike for M stay control inputs here okay but in some situations right you might not be able to explicitly solve for X dot so a more generic form here so maybe again we should write this here this is the explicit form but a merchant generic form of a set of differential equations for a general system would be something like this maybe let's say F can be f versus little F right so I'm going to try to be very careful with my handwriting at this point so little F versus big F so big F here now is a function of X dot X and u and this could equal 0 right ok so again X dot X and u these are the same things except in this case maybe terms like cosine of X 1 dot or something like that where you weren't you weren't easily able to isolate X dot and solve for that vector X dot all by itself right so this here is is another form of the of the set of ordinary differential equations right this here is sort of the implicit form of these differential equations right and again there's not a big difference here it's just a slight difference way to think about it you can easily see that if you have an explicit form it's really easy to write the implicit form in this right basically just move X dot to the other side so now what you have is if I move X dot to the other side you could write this first equation as basically 0 equals little f of X you minus X dot right so this thing here this is big capital f of X dot X and you write so if you have an explicit form of the OD es it's super duper easy to write the implicit form but again in some cases when you write your equations of motion you might end up and be stuck here at this implicit form of the ordinary differential equations right now the way we can think about this is if you think about this this function capital F right so this is f right the implicit form right what this thing is here is capital f of X dot X and you write so all this thing does is it says okay you have to tell it what is the X dot vector you want to evaluate it at you have to tell it the X vector you want to evaluate it at and you have to tell it the U vector that you'd like to evaluate it at right we already know what the dimensions of this guy's are right so we said X dot this is basically n independent variables you can kind of think about it that way right there independent variables there are inputs right there's any inputs here this is also an inputs and this here is m inputs ok and what comes out of this this zero right here is if you about this long enough this is really an N by one set of zeros right you'd actually see it right here as well right so this capital F what comes out of this thing it should be you're picking X dot X and use that are consistent with the equations of motion what comes out of this thing is basically n zeros right so the other way you can think about this is if you don't care about what comes out if you don't care that it that it satisfies this equation right you could look at capital F it's really big f is nothing more than some function where you stick in how many how many independent variables are there's n plus n plus M so this is 2 n plus M inputs and that function Maps it or spits out right n items so our n right so that's all there is to it you could just think about this big F it's nothing more than a function of 2 n plus M independent variables okay and it's a vector valued function in the sense that it spits out and that's right though the way you could explicitly see this is we could write this right you could say f of X dot X and you write this is a vector where you have maybe cap F 1 of X dot X and you big F 2 of X dot X and you all the way down to big F n of X dot excuse me X and you write where each one of these these F eyes of X dot X and you write each one of these small subscripts they're basically these are now scalar functions but there are two n plus M inputs right there are two n plus M inputs and these little F eyes they map to just a single numbers like an arm one right so it's a scalar function but there are two n plus M independent variables right okay so if we're through with that let's do one more step I think I'm gonna squeeze it in over here just so we can get everything on one page and fortunately we ran out of room actually you know what I might be able to fit it in down here to kind of keep with our flow so if you think about this uh we could just why don't I make some other vector I'm gonna call it ADA okay ADA is I'm gonna stack up X dot right so there's n of these these things we again you just think about them these are the first n independent variables or inputs to the function then I'm going to stack up X which are the next n independent variables or inputs to the function and then finally I'm gonna stack up you write something like this so again this vector n this thing is just a 2n plus m element long vector of all the inputs right so if you do that all right the beauty of this is we could now think about this function here instead of passing in you know N and then N and then M inputs I'll just pass in this vector EDA right so I could write this as simply F of EDA right so really this it's a really Tarn easy simple convenient form right it's basically saying I have a vector valued function which has to n plus M inputs right that's all we're looking at right so all we want to do and that the thing that makes it interesting obviously is capital F right this thing is potentially complicated it's potentially nonlinear and all that stuff so all I want to think about now is if if we're comfortable with this formulation I just need to find some way to linearize or to come up with some first-order linear or affine approximation of this function big f right ok so give me a second deposit camera erase the board a little bit and then let's figure out how can I go about and linearize just this function which has to n plus M inputs okay so what we're now going to do the way we're gonna linearize is we are going to live arise via the Taylor series expansion right and we actually had a video on this right where we talked about what is the Taylor series expansion and I just want to quickly recap the idea so the idea is let's say that this cat this function f is really simple it's only a function of one variable okay so it could be something complicated right it could be some nonlinear function like this right and what Taylor series said is that okay you can actually recreate this function this complicated function by picking an expansion point let's call it X naught okay and then what I need to do is you could write the Taylor series expansion of this function around this point right so I think what that said Taylor series said that okay this function f right it's basically a bunch of terms it's f at x0 right plus the derivative of F right so partial of F with respect to X evaluated at the expansion point x0 right times X minus x0 right and then plus the second-order right D F squared DX squared x0 right the second derivative of the thing at that point times X minus X naught squared and then I think we have to hit this thing with a over two factorial right plus yada yada yada right you going onto this thing for infinity right that's what the Taylor series expansion said okay so let's now expand this idea right and in the Taylor series video we said that okay this was how you did this thing for a function of a single variable now the multi the multi-dimensional version of the Taylor series expansion looked very similar so now this was the one of the case right now what if we have n D right and dimensions okay well I guess in in our case maybe we should instead of saying and we really have two n plus M dimensions right okay so in that case we had to again pick an expansion point so in our case right now with our larger system of F of a de right is equal to zero right that's what I was trying to this was what we were focusing on and what I want to do is I just want to look at the F part right so let's just look at the F part I don't need to care about the equals zero portion right now okay so we said this thing I want to expand this thing using Taylor series so Taylor series said that I can make this thing equal by basically doing a very similar situation so it's F at a de not ok so again this is now my expansion point right the place that I'd like to expand this around this plus now what do we have we have D F a not D a de right times a 2 minus a 2 not all right that was this term here plus and all this other stuff let's just call all of these things higher order terms right so plus higher order terms right that's what we had okay okay let's think about this a little bit okay so first of all let's look at Aidan not what is ADA not so in this case ADA not all right again ADA was just my function of X dot X and u right so all I'm saying is I want to pick this vector at some location which I'm calling ADA not so this is X dot naught u naught and X dot naught X naught and you not right so again all this is is 2 n plus M numbers right now again maybe this is an interesting time to to come back and say that okay right now what we're doing is we're just we're just going ahead and trying to expand this function about some arbitrary point you're free to pick this about any location you want right but that being said what we usually want to do is pick up a location of ADA naught which is consistent with F of ADA not equal to 0 right because if we do that then works spanning about a point which is consistent with the dynamics of the system right technically you don't have to write you could just say all I care about is this is expanding or creating a linear approximation of this function around an arbitrary point right it's not necessarily true that any point you pick is going to equal zero right okay so that being said we can still go ahead and do this all right okay you pick some aidan on let's just analytically go ahead and expand this guy okay about this point all right so what we can do is let's actually go ahead and maybe rewrite this a little bit okay so what I mean by that is let's go back and flip out instead of using a des let's go back and again we said hey there was nothing more than these three variables stacked on top of each other so you know what I could actually rewrite this guy okay so again the left side is still F of a de all right this is the actual thing right is equal to F okay okay instead of writing eight I'm gonna regulate this X dot naught X naught u naught right okay plus okay and now his word gets interesting this guy right here right this is a big giant Jacobian matrix right of partial derivatives of first-order partial derivatives of every function fi with respect to every a de like an A - J right so this here if I were to write this out as just a de or using X dot naught X naught and you not I'll just write it down and then I think I'm just blabbering at this point so maybe I should just write this and we could say this is basically DF the X dot naught X naught u naught D I'm going to take the derivative with respect to the X dot independent variables here this and then I'm going to multiply this now by X dot minus X dot naught okay this okay plus now I'm going to write come down here the second this is now the D F X dot X oops sorry no dot there you know what DX this right times X minus X naught okay and then finally plus D big F X dot X u and all of these are knots all right D u right times u minus u naught okay and then plus higher order terms okay so let's just take a quick look at this real fast to make sure we're all on the same page so this guy in dashed blue right won't tell you what let's look at the thing that's outside the - blue this term right here I think everyone will agree this is the exact same thing as this right I've just rewritten it again instead of compressing ADA into one giant vector I've just rewritten it down said okay ADA is really X dot X and you alright okay and now the term in blue is basically the same thing as this term in blue right I just again expanded it and rewrote it in terms of each one of them so again what this thing right here was this was a giant Jacobian matrix evaluated at the point of interest at our Taylor series expansion point of a - not right so this first term here right is that Jacobian matrix as well but I'm first I'm gonna capture only the dependence of the Jacobian on the first n independent variables which happened to be X dot so this is really here this is sensitivity with respect to variations in the parameter X dot all right and again this is a n element vector right ok then I take that element vector and I have to multiply it by the perturbation in the first n independent variables all right which is X dot minus X dot not right okay then this term right here this is now what its sensitivity with respect to the next and independent variables which happen to be X right which again is just n variables right and again well I I missed a knot here all right I got to make sure everything has a knot because this year again it's the Jacobian evaluated at the expansion point X dot naught X naught u naught which is the same thing as the expansion point of a to not right okay and then finally this turned down here this is now sensitivity with respect to the last M States which really is you right ok great ok so that's all this is this form is a little bit nicer to deal with here because we can consider each one of these three these are like sort of sub Jacobian matrices if you want to think about them this way right but this is really it's this one is measuring the sensitivity of the function the big f with respect to only variations in X dot this one measures the variations or the sensitivity of this guy with respect to X and finally this third Jacobian matrix red your sensitivities with respect to you and in fact you know what might be more helpful is less explicitly write out what do each one of these Jacobian matrices look like ok so let's do that over here and I can erase all this because I don't think we need the one B thing any longer ok so ok let's look at this first term right up here this thing D f of X dot naught X naught u naught D X dot right what is this thing right so what this is is this is a n by n matrix and let me see if I can write this like this ok the first entry right here this is going to be d f1 right remember it's the the first scalar function within the big vector valued function f this guy [Applause] partial D X whoops hold on wait look like I want to make this clean sorry D X dot one okay and again notice I don't have a bar on this right because all this this one one element what it measures right is it is capturing the sensitivity of the first function with respect to the first independent variable which another name for that is X dot one right it's not a vector it's tickets this thing is just a scalar partial derivative right okay and then the next element over here it's a similar it's now DF 1 DX dot 2 right and etc etcetera all the way down to D F 1 D X dot n ok that's what the first row looks like right so the first row of this thing is you just take the first function capital f 1 and all you do is you take its partial derivative with respect to the first n independent inputs which are X dot 1 X dot to all the way down to X dot n right then you go down here to the second row so the second row now you do the same thing for the second function so this is now DF 2 D X dot one and then DF 2 DX dot two all the way down to DF 2 D X dot 2 right and etc etc right you continue this pattern until you finally get to what you get D F n D X dot 1 DF n DX dot 2 all the way to D F n D X dot n right you take all those partial derivatives right so that's what you do and then we got to be real careful remember about this notation here of X naught or X naught X naught and you not means you take all of these partial derivatives right this thing is it's big it's also kind of ugly but you evaluate it at the end of the day once you're all done you evaluate this thing at Aidid R X dot equals 2 X dot naught at X equals to X naught and you at U naught yikes so that's kind of oh that's kind of a mouthful but that's what this ends up with right what you end up with at the end of the day is let's look at this this thing is a n by n matrix right that's what you end up with at the end of the day okay and actually let me let me scoot this up so we can see oh maybe I won't see where that's we use this thing maybe I'll write this down here this is an N by n matrix okay okay so great that is this term right here check box we got we understand what this thing looks like right it's this square matrix where it is just filled with partial derivatives with respect to the first n inputs to this function to the function big f right okay so now let's take a look at this one right here so this next one here which is we wrote here as DF X dot naught X naught u naught D X right what is this thing okay so again as you can probably imagine it looks really darn similar okay in fact I can probably just sketch out the shape of this thing right off the bat oh yeah sorry hold on yeah no that's that's fine okay okay okay so this first element right here this is going to be d f 1 DX 1 so again notice there's no dot here right there's a dot on this term but there's no dot here okay and similarly DF 1 DX 2 all the way to DF 1 DX n right and then you again follow this pattern down this is DF 2 DX 1 DF 2 DX 2 all the way down to DF 2 DX N and you follow this pattern down until you get to DF n DX 1 DF n DX 2 all the way down to DF n DX n right so you take all the partial derivatives of all of those n functions with respect to the middle and so this is now I guess it's it's maybe I'm blabbing too much I think everyone sees what we're doing here right is you take the derivatives with respect to X 1 X 2 all the way to X n right not the dots like we did above right and then at the end of the day once you take all those partial derivatives you end up with some color there's some ugly expression you evaluate this thing again at the same location right you're still evaluating this this Jacobian matrix at this particular location right so you gonna have to evaluate this thing that next time not X everywhere you see an X plug in whatever values for X naught you have and wherever you see you you plug in the values of U naught right ok so you end up again with an N by n matrix okay check box we got this thing nailed down right this this this second term so the last thing we got to do is this one again it's not complicated it's really darn similar to what we've got going on over here so this last term this matrix here right I need to see what is D F evaluated at our expansion point with respect to you right so again it is just a matrix right which captures all of these are the function sensitivities with respect to the lab em independent variables another Wharton name for that is just the vector you write so this first entry this D F 1 D u 1 and then again D F 1 D you two all the way down to now what's interesting about this is you end up with D F 1 D u M right because this element this vector u there's only m elements of this thing right this thing was a member of our M right okay so we end up with this let's go ahead and draw this guy out a little bit okay and then the next row of this thing is d F 2d u 1 DF 2 D you two all the way down to D F 2d um and then again we follow this pattern all the way down until we end up with D F n D u 1 DF n D you two all the way down to D F and D u M right and then again same thing we evaluate this sucker at X dot equals x dot naught at X is equal to X naught and u is equal to u naught right so at the end of the day we end up with something which is actually this guy is not square right those two matrices were square but this one comes out to an N by M right there are n rows because there's n functions and there's M columns because we only have M control inputs right so we end up with an N by n matrix okay okay you know what let's give these matrices names so we don't have to keep writing them all down let's call this first one here this matrix I'm going to call this thing that e matrix okay and then this guy right here which measures the sensitivity of the function how it varies with respect to the States let's call this the a prime matrix okay and then this matrix here which says in measure the sensitivity or the function with respect to the with respect to the control inputs let's call this thing the B prime matrix right okay so if we do that this actually gets a lot more simple in fact you know what we could say is you know this thing right here this matrix I think we called it e all right let me just write a big E in here this is just e and then this matrix right here I think we call this a prime all right and makes another let me erase all this my annotation so we can just get to the equation okay so I get rid of this thing okay this thing is a prime and this thing right here was what we call B Prime great okay so all right let's um let's rewrite and tell you what look look yeah let me let me rewrite this thing actually ah let me rewrite this whole thing because I think it it will make our life a little bit easier if we have a little this laid out in a little bit better fashion okay so we ended up with what we said f of ADA right this is what we are doing this is what the function and then our Taylor series expansion about this said that we had where did I have this thing whoops oh here okay this was F of ADA naught plus e times this was X dot minus X dot not right plus a prime X minus X not plus B prime u minus u naught okay great or plus higher order terms right okay okay now uh this is pretty darn handy a couple things we can note um why don't I move whoops why don't have a bar here uh sorry okay a couple things we can do here first off if we had chosen an ADA not here that obeyed our or was consistent with our with our dynamics right so if if ADA not was consistent with dynamics what do we have we had in the beginning we said that okay that meant that it satisfied this right big fat zero so this if we're expanding about a point that's consistent about the dynamics this is just big fat zero right okay if not if this thing was nonzero we could move this to the other side so we could write this as f of ADA all right you actually you know what let's leave it like this this is maybe the easier way to go about this okay okay so this is zero okay what is this left side right in fact remember we said the the left side don't we have that F of ADA is equal to zero right that was our rig this is our original implicit set of equations right this was our original nonlinear differential equation so this is also zero right so we could rewrite this whole thing as zero equals E times and now let's look at this term right here let's call this if you think about this long enough this is basically saying how far away are these states or are these independent variables from my trim point of interest so this is almost like a delta or a perturbation of how far away did you go from the point X dot naught I'm gonna call this let's call this thing how about Delta X dot right because it's talking about how far away did you go from X dot naught right okay similarly this term right here this is how far away did you perturb er how far did you move from the the expansion point X naught let's call that a Delta or perturbation again this is like Delta X right and finally this term right here this is Delta U right okay and yeah that's that's fine okay so this is e times Delta X dot right plus a prime Delta X plus B prime times Delta u plus higher order terms in the Taylor series expansion right so again right now we've actually I can keep writing an equal sign because we haven't made any approximations of course this higher order term thing goes to infinity so this is not terribly helpful so what I think we should do now is let's make an approximation and drop all these so let's drop these higher order terms and assume there's zero okay if that's the case I now have to write an approximate sign alright but now we have approximate approximately equal to e times Delta X dot plus a prime Delta X plus B prime Delta u okay let's go ahead and solve for Delta X dot right so move this thing to the other side and we have what is this yeah did it yep must be negative yeah so actually tell you what let's look let's move the other junk to the other side so we can write this lovely I think I'm rambling at this point so e times Delta X dot is equal to negative a prime Delta X plus all right minus B prime Delta u right okay now to isolate this Delta X dot I just need to let multiply everything by the inverse of e ax right okay so if I do that he hikes I'm almost out of space but what gosh we're so close I tell you let's come over here then let me see again I apologize for the whiteboard usage space but okay isolating that Delta X dot is now equal to negative e inverse shucks I'm running out yeah okay sorry I tell you what let me let me do this down here uh-huh I kind of want to keep this here this is bugging me we're literally out of space okay let me see if I can I can cram it in here okay so we have Delta X dot is equal to negative e inverse times a prime Delta X minus e inverse B prime Delta u right why don't I call this quantity here this thing this negative e inverse a prime call that a right because this is an N by n at times an N by n so you get some other N by n right similarly let's call this term here this negative e inverse B prime let's call that another matrix B okay so if you do that we end up with Delta X dot is equal to a delta X plus B Delta u right and what BAM look at this this is our linear system that we're looking for where you're a matrix right so the a matrix is just going to be negative e inverse a prime and the B matrix is going to be negative e inverse B prime right this is a and again I apologize that we crammed it in this little spot because this is that this is what we're looking for right if you look at this this is your familiar linear time-invariant system right it says Delta X dot is equal to a delta X plus bu or if you want to drop this Delta notation it just says X dot is equal to ax plus bu this is our linear system and what's super important about this is here are the formulas of how to get it and we started with an implicit nonlin your set of ordinary differential equations what we can then do is you go through these Jacobian matrix calculations right here's how you get e here's how you get a prime here's how you get B prime then you just do a bunch of inverses and you multiply the things together and you basically will get the a matrix and the B matrix here so that's the theory behind the linearization okay the next question is we saw that all this hinges down to is it's nothing more than calculating the e a prime and B matrices boiled out one step further that is nothing more than calculating all of these partial derivatives there's a lot of partial derivatives here right this is an N by n matrix so I think there's what there's n squared terms here plus another N squared terms here plus n times M terms so there's a lot of partial derivatives we're gonna have to calculate okay so that's what I want to talk about Nexis how do we calculate all of these partial derivatives right there's many different ways to go about that but as soon as we're able to calculate all of these partial derivatives we should be in business okay so give me a second to pause the camera erase the board and we'll be back in just a second all right so I reached most of the board but I kept up the description of the a prime matrix just so we have something to reference so what I want to talk about now is okay how are we actually going to go about calculating the e a prime and B prime matrix right so we need to figure out how to calculate these Jacobian matrices now the thing that you can probably see right off the bat like we discussed is okay your first impression might be to say okay these are just a whole bunch of partial derivatives so all I need to do is just go and take a whole bunch of partial derivatives right that's totally fine and that's that's the analytical approach right you can just say if you had an analytical description of all of the functions F 1 F 2 all the way down to F N and they were simple enough you could just go over to your favorite symbolic manipulation toolkit or heck calculate my head if you want but it's just a bunch of partial derivatives right now what I do want to talk about is what happens though if this is not available to you so what if analytically calculating these partials is not an option so for example the function that describes the the dynamic fear system what if it's a very complicated function or what if it's a black box like we talked about where you don't have an analytical description I can't fit this function to cosine science exponents polynomials things that have analytical derivatives right so in that case we need to start thinking about how can I get these numerically right so we may come up with some numerical approaches to go ahead and calculate these functions okay so the first way I want to look at this is we are going to do this by numerically calculating each one of these individual partial derivatives using a basic finite difference method in fact we had a dedicated discussion and video on how to numerically calculate partial derivative so that's this video over here so make sure now if you haven't seen this video on numerically calculating partial derivatives pause the current video right now go watch that video make sure you understand it because I'm going to leverage that idea that we developed in that video right here okay so in that video what we did is we figured out how to numerically calculate the partial derivatives of a scalar function of potentially n independent variables right if you look at this right now and say if we just isolate a single row of one of these okay so in a circle this may be in blue okay just focus on this row of say VA prime matrix right if you look at this long enough what this is doing right this is really this row is calculating or asking what is the partial derivative of what is the sensitivity of this scalar function f1 so Row 1 right is basically looking at it focuses on only f1 of X dot X and you right so this function after one right so if you recall we said that F 1 this is just a function that has two n plus M inputs and it maps to a scalar output right this is exactly the situation we looked at it when we were talking about how to calculate the partial derivatives of this scalar function with respect to all of its inputs right that's basically keep talking about numerically calculating the gradient of this function and we saw we could do that using the symmetric difference quotient approach so that's all we need to do right here for both the e a prime and B matrix right we basically are going to go ahead and and do that so if you go ahead and basically use the symmetric ah keep somebody's wrong submit trick difference quotient right what we end up with is we end up with an expression of each one of these entries so for example let's look at in this case the a prime matrix so I'm going to write that down here so so a prime the I J entry here of this guy is what is this thing so this is basically it's DF I right it's the I so the yeah the Row is the il iment right so it's this thing D X J right that's this thing and this is here evaluated at the expansion point of interest 8 and not right that's what this element of this matrix captures right so if that's what we're trying to compute we can basically now say this is approximately equal to that rise over run formula that and basically apply the symmetric difference quotient right here to try to get an expression for this specific partial derivative right so if you do that we see that this is basically going to be the rise is going to be the function f I write the I function evaluated ad this is where we got to be a little careful right so now I'm gonna evaluate this guy at X dot not right and now I'm gonna perturb X here and we're gonna perturb it in the j DH element in the positive direction okay and then you not this okay - f I oh gosh I ran out of space a I should have should have made this a little bit easier maybe maybe we can fit this all in X dot not X I J - and then you not yikes barely fits okay okay so that's the rise and then over the run here is two times whatever the perturbation distance is so here Delta X I J okay okay so in this case again with the a prime matrix all we're doing is we are only perturbing the appropriate independent variable so the appropriate independent variable right a prime is measuring the functions sensitivity to the X to the States right X here the not not X dot right just the states X right so in this case we only perturb if I want the J element right the J column of this matrix you only perturb the Jade state so this thing right here this X let's write it over here where X I j+ what this thing looks like right this is x1 not x2 not all the way down to X J not but then you perturb XJ in the positive Delta X IJ direction or in this direction by this perturbation amount right and this all the way goes down to X and not okay with this and then I'm gonna transpose it right because this because this vector should be a vertical vector right the the we kind of interpret it right so again all you're doing is you're carefully your your basic this vector is basically the expansion point all you're doing is you're just perturbing only the J the state by the distance Delta X IJ and this is some small number right okay and then this other term here the the part of the symmetric difference so this is the perturbation in the forward direction I perturb the Jade state forward by Delta X I this term right here this this X IJ minus this is in the negative direction right so here X I J minus this is again it's very very darn similar x1 not x2 not all the way down to X J not but now I subtract off Delta X IJ and all the way down to X and not and then transpose it there we go okay so that's the a prime matrix the e prime entries are very similar so let's write that up here so here e sorry it's not a prime is just e e IJ the J faith ro j DH column of this here right what that is capturing again it should be the function it's the ithe function right that's the row again this things dependence on X dot J all right and again we evaluate this sucker at eight and not that's what that element of the matrix is so trying to use the symmetric difference quotient is basically saying I need to now do a very similar idea right it's I need to evaluate the I function at X dot i J plus write and then X naught u naught this that's a forward direction and then minus fi X IJ minus and then X naught u naught okay there we go that's the rise all over the run of again two Delta X dot I J write it again this perturbation distance can be different than that right and again the this perturbation vector this guy and this guy okay let's just write them down so here so we're X dot IJ plus is X dot oops sorry X dot one not X dot two not all the way down to X J not dot all right and then I perturb in the Delta X I J dot yeah I'm sorry this is getting a little bit messy here but I think okay maybe let me erase these these arrows just so we leave only the things we care about you got to transpose that okay so that's a perturbation vector so again this is this perturbation vector where all I'm doing is I'm perturbing in the j DH independent variable right so all i have to do is create this this vector is basically it's basically this thing X dot not but all I'm doing is I'm perturbing the j DH element forward by Delta X dot IJ right and then similarly the minus term right so X I J dot minus is basically the same thing x1 not this should have dots x2 not dot da-da-da-da-da all the way to X dot J naught minus Delta X Phi J dot da-da-da-da-da X and not dot likes yeah transpose okay sorry that got a little more crammed that I wanted to okay but there we go here's the e matrix here's the a prime matrix so the B prime matrix is again very very darn similar we are just going to apply the symmetric difference quotient to this guy so again the IJ element of this guy is going to look like again what it should be right analytically or theoretically this is the sensitivity of the eighth function with respect to the jet control input right again that is evaluated about our expansion point that's what this is and then we apply the symmetric it's quotient okay so we need to now perturb the function the eighth function okay in the only in the controls now so this so these are the expansion points but now I need to perturb the control so I take you I J plus write a perturbed this in the positive direction and then minus the function evaluated at again the expansion point almost and then here's you I J - okay so there we go there's our rise over our run is a very similar expression 2 times Delta u I J ok great and then just for completeness let's write down exactly what is the perturbation vectors here it's gonna look very dark very very similar so here you I J plus alright it's basically you one not you two not all the way down to you J naught and you perturb this sucker in the positive Delta u IJ direction and this now goes only up to you M right there's only M controls okay this thing transpose okay huh almost done lastly we got you IJ - all right so this is you one not you - not all the way down to you J not now - Delta UI J all the way down to UM naught transpose yes there we go ha so this is basically the symmetric difference quotient so basically if you think about this long it's just a few for loops this email just have to go I maybe we should maybe we should have written that what are the ranges of I and J so in this case obviously I and J are our integers but I and J go from so I goes from 1 to n right and J goes from 1 to N right because I going from 1 to N means there are I function alright sorry n right here's one - well I guess I'm pointing at the wrong thing but I think you get what I'm saying there's n functions and J is how many of the independent variables are we perturbing in this case so there's n of them we're doing here now for the a IJ this expression here I also runs from 1 to N and so does JJ also runs from 1 to N ok but then the B matrix we get a slightly different expression right so I runs from 1 to n right there's still n equations so there's still F 1 F 2 F 3 all the way down to F n but there's only M control so in this case J runs from 1 to M ok so this one you would stop this for loop a little bit earlier right but long story short this is how this is how you would get it um and always the other thing that I maybe want to point out before we think about trying to how to implement this is we talked about this during the lecture on numerically calculating partial derivatives right but the perturbation distances that you choose they could be different depending on the function and how it depends on each one of these variables right if this if the function has no dependence or a very slow dependence on one of these independent variables you can choose its corresponding perturbation distance to be to be large right but if it changes fast you have to make that smaller so that's why we're using this notation for how all of these perturbation distances can be can be individualized right so you can make this thing that's why I'm using two indices here right so you can use a perturbation distance which is associated with the i'f function and the Jade state all right so this is a it's sort of the most general Terfel format and formulation but in practice right just make all of these numbers the exact same small numbers probably a good way to start I mean I wouldn't use machine precision but you know 1 e to the minus 8 or something like that might be good so maybe we'll just make a quick note here I'll just write it down just for completeness no in practice you know using Delta X dot IJ just make that thing the same number as Delta X IJ which should be the same number as Delta u IJ right just make them some small number and you all let you play with that what is a good small number to get to get reasonable results right okay so that's the game plan at least that's in theory right we haven't applied this to any particular model this is just a theory of saying okay we need the e the a prime and the B prime matrices all those things are there just Jacobian matrices which are just filled with partial derivatives all we're doing with this scheme here is we're just we're brute forcing it in the sense of we're just going through and iteratively calculating the partial derivative of every single element in this matrix using the symmetric difference quotient approach okay so let's jump over to a system now give me a second to erase the board and then how can we apply this to a complicated dynamic system so the yeah let me erase the board then we'll come back and look at an example alright so the specific system that I want to look at linearizing using this symmetric difference quotient numerical method here is a nonlinear six degree of freedom rigid body aircraft model in particular I want to use this research civil aircraft model or the Arkham model again this is maybe a time to pause if you don't know what I'm talking about when I say the our camp model please watch these videos these are videos where we talked about developing the equations of motion for this aircraft and we saw that they were very complicated but that being said we were able to build a MATLAB Simulink implementation of this model and more importantly we were able to derive and build a MATLAB function which was able to get all the equations of motion in explicit form and basically return the state derivative if given the state in the control so again all of that information is captured in these two videos so make sure you've watched those two videos and you have a good understanding of what this model does so that we can go ahead and verge it okay so that being said the aircraft model that we developed was a couple of things so the thing that we made in that lecture was we made this thing called our cam model dot M right there was this mallet function where if you gave it the states and the controls right this thing would spit out X dot right this is that that function that we spend a good bit of time developing okay so that was all fine and dandy but we did say that this was here in explicit form right and you can see it's explicit because that's what this function did write it explicitly solve for X dot so in order to use the meant the the more general method that we just talked about about being able to linear as an implicit set of differential equations it's really easy because we already have this right so this guy is basically the explicit form of the Archaea model so all we need to do is come up with another function so I'm gonna use MATLAB speak right now to define this function so we just need to say something like function maybe I'm going to say F Val for the function value and this is just going to be I'm gonna call this something different like our cam model implicit okay and then this thing you're gonna give it X dot I'm gonna give it X and you okay and that's all this function the function I can almost make this in one line right the implicit version of the function is basically the function value of this X implicit version it's basically what it's basically you just say our kam underscore model you give it X and you write and then minus X dot alright that's semicolon that's literally it so the function it's it's a one-liner you make this M file which is literally well I guess it's technically two lines but right then is basically implementing capital f of X dot X and you write that's what this function implements right and the way it does this is it goes it says okay this is basically lowercase F of X you write minus X dot kind of the same thing we did at the very beginning of the of the discussion tonight right so this is basically the call to the our cam model right here right that's this and then all we're doing is subtracting off X dot right so that's the first thing we're gonna have to do is if you build this this M file right this is basically encompassing the entire complicated dynamics and equations of motion for this six degree of freedom rigid body nonlinear aircraft model so that's this big capital F right we basically have made a matlab implementation of this big capital f function by leveraging all the work we already did all of the hard dynamics were captured in this thing this and this our cam model M right which described the the vehicle dynamics okay so really the modification to make it in a format that's usable for our implicit numerical some numerical differentiation via the symmetric difference quotient approach it's it's really pretty darn simple so we'll go over to the board and we will go over and do this in just a second ok so we're almost there right so if you remember the goal of the of the numerical integration scheme is to develop and get right so our goal right now is to numerically calculate e a prime and B prime right and again just refresh your memory right so all this the e matrix was the matrix was basically partial of big f with respect to the first n independent inputs to this function right but this thing was evaluated at X whoopsies maybe what we should do here is this is evaluated at a de is equal to a to not right and same thing with the a-prime and b-prime maybe let me give myself a little more room so I don't crowd myself right so the a prime matrix was basically D capital F D X all right and then this thing is evaluated at a de is equal to a and a naught and finally B Prime right is D big f D you evaluate it at a de is equal to a de not right so all we're gonna have to do with our with our scheme is we're gonna have to differentiate this function we're gonna have to do a numeric differentiation to get the e matrix right we are going to perturb these first and States or the first X dot independent variables are independent inputs to the function right that gives us e Prime so here are sorry this is e right and that the a prime is perturbing this guy right and then the B prime is perturbing this set of inputs ok the last bit that we have to do is we have to do this at a specific a to not right we have to pick an expansion point that we want to linearize about ok so the specific point that i want to use again we talked about this way at the beginning you want to choose a point that is consistent with the dynamics of the system right so I want to choose I want to put the vehicle or the the system at some value of X dot X and u which corresponds to some point of interest right so the point of interest that I want is this steady state straight and level flight i want to linearize my airplane as its flying along straight and level going around at a specific speed now again we had a dedicated lecture talking about how to find the x and the you and actually the X dot as well how to find this value of a de where the plane is in this configuration so again that's this lecture over here so if you haven't watched this video again take a moment pause this current discussion make sure you watch it view that video where we talked about how to trim the meet vehicle and find specifically the X U and the X dot which corresponds to straight and level flight okay so I'll assume you did that and tell you what looks let's just write it down so we're all on the same page so the trim X maybe let's write this well we can start with the trim X dot right let's let's do that how about X dot not what is the trim value for X dot and that is pretty easy right because this condition of steady state right we saw in that earlier video that in order to have steady state you got to have X dot is equal to zero okay so the aircraft is the states are not changing so this is a big fat vector of nine zero so here one two three four six seven eight nine okay there we go that's X naught we need now X I'm sorry said that's X dot not again there are subscripts are live getting out of control but I think we're almost there okay now what is X not okay X naught is the state's at steady state straight and level flight so I think in our case this was okay what was you it was eighty four point nine nine oh five and then zero one point two seven one three and then P Q and R over L zero and so was Phi was zero but theta was non zero I think this was zero point zero one five over radians and then zero okay there we go okay so here was the state to get you trimmed straight level flight at around 85 m/s airspeed right so here's UVW P QR v theta sigh okay and then finally the controls that would get you that condition that's you not this was what was it it was 0 and then minus zero point one seven eight oh and then zero and then the two throttles were matched here in 0.0821 0.0821 okay so there we go here's eight and not all right it's these nine plus nine plus five values right so this is the trim point that I want to linearize about we now have the the model here and now we're gonna combine this with what we just discussed and we're gonna go and let's go ahead and jump over to my lab and see how can we write ourselves a function right that is going to go about and calculate the e a prime B prime matrices by doing the numerical the numerical approximation of the partial derivatives of each one of those matrices using the symmetric difference quotient and we want to do it around this this point right around steady-state straight and level flight so okay um that's the game plan let's jump over to MATLAB and see how we're gonna write this this function that will do this linearization alright so because we've done all of this work already on the actual implementation isn't so bad so again here's my our camp model M this is the large function that we built earlier which basically contains all the dynamics for this aircraft right so this thing is again what we spent a lot of time on previously and this embeds all of the complicated dynamics associated with this particular system so then our our cam model implicit function like we just talked about this is basically got two lines of code right in this case now our can model the implicit version right you give it X dot X and U and we do just like we did on the board right all you do is you take here's lowercase F of X you minus X dot right and this is the implicit version of those exact same dynamics now here's where it gets interesting I've made myself a function that I'm calling implicit Lin mod this is what is going to do all of that calculating the partials of every single of the e prime the sorry the e the a prime of the b matrix is going to calculate this using that symmetric difference quotient technique we just talked about so what you do is you pass this function a function handle so we're going to pass it a function handle to the Archaea model implicit and say that that's the function that I want you to interrogate and calculate the the matrices for me okay so to do that you have to pass it the point that you're interested in again so here we go here's the effectively these three inputs are eight and not right this is the X dot not the X not in the unit so these are the basically the expansion point for my Taylor series and then what I'm going to ask for is like we talked about on the board is you need to choose and define what are all of the perturbation the perturbation distances for every single variable in every single function so these guys are also matrices of the appropriate size and all of them you know to make again to make your life easy every these could be matrices of the appropriate size but they could all be filled with the same small number heck I mean if you want to you could just make instead of passing in um three separate matrices like I've shown here you could just pass in a single specific number which is that small perturbation distance and you use that same perturbation distance for every single function in every single variable but again I've tried to make this general so if you if I wanted to perturb different functions in different directions or sorry different functions different distances you could go ahead and do that here using this framework so again that's the function prototype right you pass in all that information and now what it does is it goes through and rips through and computes the e the a prime and the B prime matrices okay so let's just take a quick gander at this what it looks like so the first thing I do is I figure out how many states there are right you can take the length of either X dot naught or X naught both of those should be the same thing and then I'm gonna find out the number of controls by taking the length of U naught okay then all I have to do is a start going through those matrix calculations and we saw that what we're gonna do is I'm just gonna brute force my way through it so let's start with the e matrix whoops sorry like here so we know that e matrix is going to be an N by n matrix so I'm just gonna start here and initialize the container of all zeros and then once I have this container set up all I need to do I'm gonna write myself a double for loop so I'm gonna loop through every single function in the big F matrix and I'm gonna rip through every single dependent variable well not every single I'm gonna rip through all of the dependent variables of X dot here right okay so that's what the e prime matrix is right it calculates the functions dependence on X perturbations in X dot okay so the first thing I do inside this double for loop is I extract the appropriate perturbation magnitude that I want to use for function I and independent variable J right so again I just pull it out of that that matrix that I pass in again think about this as just this is just gonna be a small number right and then all I need to do is I need to make the X dot plus and the X dot minus or I think we called it on the board it might have been X dot IJ plus and X dot IJ minus but long story short this is the perturbation vector in the positive direction so what I do is again I start and we say that okay you're gonna start at the expansion point right and then what I do is I tap on I only perturb the J's element of that guy and I perturb it by the perturbation distance in the positive direction that's why it's plus X dot right so again I'll let you look at this code and I think if you look at it long enough you'll convince yourself that what at the end of line 22 after you execute line 22 this thing that I'm calling X dot plus looks very very close to the to the expansion point X dot ah not accept its perturbed only in the appropriate variable by the appropriate direction right and then I do the exact same thing in the minus direction right I start and say okay let's start at the perturbation at the expansion point right and then I'm only gonna perturb the the row of that matrix or the row of that vector which corresponds to the appropriate independent variable that I'm trying to pretend now I'm going to perturb this in the minus direction right so this is most of the work once I do that we can go ahead and evaluate the function at the positive two location now what I'm doing here is the line 26 here is I am going to call my function so in effect in fact I'm calling our cam model implicit here right so I'm gonna call that by passing in this perturbation vector and notice here I'm only perturbing again for the e matrix all I care about is perturbing the X dot those those set of independent variables that are associated with X dot so I leave X naught and u naught at the expansion point okay so this thing this function R can't model implicit it's actually gonna give me an entire function so this F is going to be an N by 1 vector ok now in this thing um then my double for loop I'm going through um let's think kind of ways how I'm just I should say this in this in this double for loop I'm the outer loop is is the function value so what I'm doing here is really this is actually slightly inefficient code but I wrote it this way because I want to be able to easily understand what we're doing what I'm trying to get at is line 26 if you think about this long enough you're actually perturbing and you're getting the capital F perturbed in the positive direction associated with variable J but you're getting function f1 f2 f3 all the way down to f9 so really this is actually doing a fair bit of work here that we could leverage if you wanted to optimize this code but what I'm doing again to keep this thing understandable and make it fit what we did on the board I am actually going to only keep one element of this so I'm throwing away eight bits of information and I'm only keeping one of them either I'm only keeping the one we're focused on right now okay because I'm only trying to calculate row I of the e matrix really you could calculate multiple of them simultaneously again if you wanted to optimize this code but I think I'm belaboring the point at this point so I just want to say that that's what I'm doing right here I'm throwing away information that we could have used for the sake of understandability so this thing right here at 27 this is just a scalar value right this is basically F at X dot plus X naught and you not okay and same thing I'm doing the same thing for the - uh erection right so the rise is basically this - this right and you can see me calculating that right here here's the rise and then I divide it by the run so right here line 34 all I'm doing is I'm calculating that partial at location at row I column J and I'm calculating it via that symmetric difference quotient approximation right so you see that by the time you blast through this double for loop you have filled up the entire email takes right now let's keep scrolling down and take a look at this the a prime matrix it's the exact same thing almost right so again initializing container the only difference right is now instead of perturbing X dot right you're now perturbing X so now I do the exact same thing all right you start here at at the expansion point and then you perturb in the positive direction and you perturb in the in the minus direction right so it's virtually identical so I do the exact same thing calculate it and store it in the matrix again so there's nothing to really worry about here it's the exact same thing then just going down to the B prime B prime again is also very similar again you just kind of notice the only difference right now is the number of columns of B only goes up to M it doesn't go all the way up to N it goes to M Mike M Mike right not n November okay so again exact same thing so there's no need to really belabor the point I think you see how this works so this function you can see it's just a series of three sets of double for loops which are gonna go through and calculate those three matrices okay so now what we can do is we're basically ready to rock so here's my script good old friends to clear CLC close all and then all I'm gonna do is I'm gonna define the expansion point that we're interested in here so again here's the nine zeroes for X dot naught here is the state to be in straight and level steady state flight and here's the the controls to be in trimmed steady state straight and level flight so then what I'm gonna do right here in line 45 through 47 you can see I'm basically you know for all this work of setting up all those different perturbation distances I again use the exact same value for every single perturbation distance for all functions and for all variables I'm using you know 10 e - 12 or 8 some small number and then all I do is I call my function that we generated right this implicit Lin mod and I say go ahead and linearize the the implicit set of equations about the expansion point of steady state straight and level flight and basically use 10 to the minus 12 is a perturbation distances for all of them and give me the a or the e a prime B prime matrices once I get that I'm going to go ahead and calculate a the a matrix as we talked about on the board right it's just negative e inverse AP and then similarly blah blah blah all that kind of good stuff then I'm just gonna print them out to the screen and then we're ready to rock okay so if I let this guy go let's go ahead and run this thing right now look at that it just blasted through and it had no problem let's go ahead and look at the results um here we go well actually here to tell you what let me let me move this over and I'll run this again maybe we can get a to show up on one row let's run there we go okay so here's my a matrix right it's a it's a it's a nine by nine a matrix and here's my B matrix it's a nine by five so this is my linear system representation of the aircraft um maybe before we leave this what we should do is let me print out a as well and if you think about this look at this e is basically the negative identity matrix and that makes sense right if you think about how we set up the equations of motion because our equations of motion were able to be written in explicit form you can see here's basically the negative identity matrix coming from right here all right so this is pretty darn awesome we basically have gone ahead and we have linearized the system this a and B matrix is my linear system of my aircraft about steady state straight and level flight at the condition we specified so ah before we dig too far into this let's think about this let's let's maybe look at this matrix and ask if it makes sense right one thing to notice is look at this in the a matrix this last column is all zeroes right what that means is it's basically telling us that the function has no dependence on whatever the ninth state is and if you remember a state vector right it was uvw PQ R Phi say the Phi Theta sy sy being the heading angle of the aircraft so this is interesting it's basically telling us that no matter what heading angle the aircraft is at it behaves the same and if you think about that that makes perfectly good sense right the aircraft should behave the same if it's flying north or if it's flying east or if it's flying in any different direction it should have the exact same dynamics right okay so now that we've got these a and B matrices maybe what we should do is there's more than one way to go about calculating and doing this numerical linearization so tell you what let's let's pause the camera and let's talk about another technique to go about obtaining the similar results all right so that was pretty exciting we saw that that method worked in the sense that we definitely got an A and a B matrix but now we want a validate then make sure that we can get that similar results using different techniques just to make sure that we're confident with the approach so if you remember there's more than one way to linearize a nonlinear dynamic system particularly if you already have this in matlab and simulink in fact this video right is one of our previous discussions where we showed how you can use two techniques in MATLAB Simulink right you can either use the linear analysis tool to go ahead and try to perform a very similar perturbation analysis but have Simulink do all the heavy lifting and the work for you so we're going to go explore that and then we also discuss in that exact same video also how to use MATLAB Lin mod function which basically calculates a linear model from a set from a nonlinear system ok so if you haven't seen this video again you're gonna have to pause watch that video so you can get all the content and understand the background and how to use this because what I want to do now is I want to basically use these techniques so we're gonna go really quick and we're just going to apply these ideas to our particular system just to double-check ourselves to make sure that we all get the same a in the same beam matrices using any one of these three techniques right the first being the one that we just did that's our hand roll technique of using numerical approximations of partial differ the symmetric difference quotient approach then method two we're going to use the linear analysis tool and then method three is we're gonna use Lin mod let's check that of those okay so let's go try method two and method three right now in MATLAB Simulink alright so we're back I'm going to assume that you've already watched our previous discussion on how to linearize the system using the linear analysis tool and Lin mod so let's jump right into it so I've gone ahead and made a Simulink model that is appropriately set up and instrumented for analysis using the linear analysis tool again if you're unfamiliar with what is going on here feel free to check out that other video where we discuss this in detail so I'm gonna come here to kin analysis control design here we go linear analysis this is going to boot up the linear analysis tool so we've got that over here and now what I can do is let's just linearize this thing about that operating point that we already have loaded in and set up and yeah look at that it successfully got something here and made a linear system so I'm gonna drag this from the linear analysis workspace back to the MATLAB workspace and oh yeah here it is Lin sis 1 so let's check it out this is what the linear analysis tool came up with and if you look at this here here's my a matrix this looks really familiar and here's the B matrix again that looks really familiar you know we should probably do is let's look at this thing Lin sis 1 dot a and let's look at how does it compare with the a matrix we just calculated using the symmetric difference quotient approximate right so let's just look at the difference between the two right and maybe what you do is let's take the absolute value of the difference and then we'll look at the maximum difference across all elements there and this is this should be the maximum difference in any element of the a matrix and look at that 1.5 e to the negative 4 that looks pretty darn good so I would say those two two techniques came up with a very similar a matrix how about a B matrix let's let's check B there you go look at that that's an even better match so yeah this is looking really darn awesome um let's go ahead and also how about verify with Lin mod right so again I'm gonna assume that you've watched the video that we referenced and you know how to use Lin mod so I'm just gonna go ahead and type in the command right here so I'm gonna linearize the model that i have already set up that's right here let me drag it over and show you so again here's the exact same system and it is set up and instrumented and ready to be linearized using Lin mod so when I call Lin mod and I say I want to linearize around straight and level like that right it's gonna give me an a B C D matrix which I'm gonna call ALM for Lin mod blmc LM dl m so again let's run this and yep here we go we got an ALM and again that looks very similar so let's do a very similar check let's go ahead and compare the maximum difference between a Lin mod and our approach and hey there you go that's that same number 1.5 e to the negative 4 so this is looking really good let's check the whoops gotta check the B matrix there you go this is looking awesome so at this point I think we've got ourselves a pretty good linear representation of our system right so here again is our a and our B matrix so this is I'm pretty confident that this is the linear system about steady-state straight and level flight so now we're in a real good position to start analyzing that the aircraft as it flies are analyzing this system that we've linearized around this point now that we have the linear model of the system all kinds of tools come to bear right we can start talking about you know for example eigenvalue analysis so i can start asking is the airplane stable as the system stable at this location and all the other kind of fun things we can start designing linear controllers all that good stuff so that is probably the topic though for another video in fact our very next video will start talking about manipulating this linear model for this particular aircraft now that we're able to compute it alright so that's pretty exciting we saw that we were able to use three different numerical techniques to get a linear a and a linear b matrix for our system now that being said we saw they all agreed well but we were real kind of comparing numerical results to numerical results so I guess I'm not terribly surprised that they all came out agreeing with one another it would be really nice if we could compare the numerical results we just calculated with say some analytical results right if we were trying to calculate those partials and those Jacobian matrices analytically now we said earlier that the whole reason we went the numerical route was because calculating those those jacobians analytically were come what somewhat intractable because the equations were so darn ugly well that's not actually the inque the case for all the equations namely what I mean by that is let's write down the equations that we were dealing with right we said that I was trying to look at this equation right f of X dot X and you write and all we were doing was we were going ahead and trying to get ourselves a linear approximation of this function about some expansion point right now we said that for the Archaea model this was just lowercase F of X you write minus X dot right okay that was the equation now this lowercase F X U is what we developed in the Arkham model dot M right in those previous two lectures we showed that this thing is kind of ugly right it's a nine element vector valued function but really all the ugliness lived in the first six equations right because the first six equations gave us expressions for UV w PQ r dot that was where all the things like the aerodynamic forces and moments and all that kind of jazz was buried the last three were actually the Euler kinematical equations which were actually not so bad so what I mean by that is let's write this out so the left side of this okay is big F one big f 2 right all the way down to big F seven big f eight big f 9 right this is equal to over here we had little F one little F 2 all the way down to little F seven little F eight little F 9 right minus X dot 1 X dot 2 all the way down to X dot 7 X dot 8 and X dot 9 okay and again we said the ugliness lived in these top equations here this is the portion these first six equations were the ugliness right all the aerodynamic forces and moments lived in these expressions these expressions were ugly in fact we flashed them up on the the board earlier when we were talking about the Archaea model specifically we solved these were ugly now these bottom three are actually not that bad right because if we were to just examine the bottom three alright let's just go ahead and draw a line here and let's isolate these bottom ones so I want to look at just a smaller 3 by 1 element of f a big f7 big f8 big f 9 alright so this is big f 7 figure a big f 9 right and that was just our Euler kinematical equation right H of Phi times Omega B with respect to e express in the body frame that's this portion here and then I need to subtract off X dot seven x dot eight X dot nine okay all right so let's go one more step let's write out what this H matrix was and what this Omega B with respect to e is and I think that'll give us a good example and we could see what is going on so again the left side is big f7 big f8 big f9 all right it is equal to okay that H matrix it's a three by three right well maybe I'll try to just kind of delineate it like this okay so it's a 1 sine of X 7 and then times the tangent of X 8 and then this was cosine of X 7 plus tangent of X 8 and then we had a 0 cosine of X 7 negative sine of X 7 and then a 0 sine X 7 all over cosine X 8 and then this was finally cosign X 7 all over cosign X 8 all right okay so this was the H matrix and then the Omega B with respect to e right this is just x P Q and R which the same thing as X 4 X 5 X 6 right okay and then finally we need to subtract off X dot 7 X dot 8 X dot 9 right okay so look at this this actually isn't so bad so for example let's let's look at the equation for just F 8 how about okay so you can see I can multiply this through and actually this turns out to not be so ugly so f 8 here is basically just going to be it here it's this line so I guess this year was this okay so we basically have cosine of X 7 all right times X 5 then minus sine of X 7 times X 6 all right minus X dot 8 right so here look at this this is not bad at all all right this is a pretty reasonable analytical expression so for example going back to our numerical approximations right for example we said that the e let's calculate the e 8 how about 1 entry right but he ate one entry was partial of F with respect to X dot 1 right evaluated at some location a 2 not right so you can see all you got to do is take the partial of this function with respect to X dot 1 look at this I don't see X dot 1 anywhere in here so this thing is 0 right in fact you can see all of these en treats are going to be 0 until you come here to X 8 right so we see that really e 8 8 is the only entry which is nonzero and basically this is DF 8 DX dot 8 right evaluate it at a de not right and you see from this you basically get a minus 1 right so great this is easy to analytically check the email takes how about the a prime matrix so the a prime matrix right let's look at a prime again let's look at the eight one element alright of the of the a prime major because this should be d f8 D x1 right evaluate it at eight and not all right so you can see all you got to do take the derivative of this thing with respect to x1 again you see x1 doesn't show up anywhere so this thing is a big fat zero right so you can see that actually the a matrix the a prime matrix the only place that this is going to have a nonzero entry is where is in which columns the the fifth column the sixth column and the seventh column right every other entry should be zero and in fact that's what we saw when we computed our well yeah the a matrix are basically eight minus e inverse times a prime but you basically get what I'm saying right that a lot of these entries turn out to be zero the only time they're nonzero is if the function has a dependence on one of the variables you're taking the the derivative with respect to so in some cases it's actually not that bad to go ahead and try to analytically it verify so in fact the bottom three rows of the a matrix and the bottom three rows of the B matrix should actually be pretty easy to analytically verify all you need to do is write these equations out and take a whole bunch of partials right and actually costs just looking at this and you know in fact you can you could pretty much tell that the B prime matrix rows 7 to 9 I'm going to use MATLAB speak here writes of rows at the bottom three rows all columns right if you think about this long enough what is this going to be it's the bottom three rows so it's the functions dependence on you on all users right on all the use this thing is going to be basically it's a big fat set of zeros and the dimensions it should be three rows by however many controls are were in this problems I think there was five right because you see in this equation up here I don't see you showing up anywhere there not a you want to you - there's no use at all anywhere so these are going to these bottom three rows are going to be a big fat zero so this is actually pretty handy it's a good way to go about checking to make sure that the a and the B matrix align with our our analytical analysis here all right so I want to summarize what we've talked about tonight so what we've done is we've basically developed a system to numerically linearize a nonlinear dynamic system about any operating point that we choose right and this is incredibly powerful because we saw there's three different numerical ways to do this you can either use the symmetric difference quotient technique you could use the MATLAB Simulink the linear analysis tool or you could use MATLAB Simulink lin model functions and all three of them seem to give you very similar results okay now that we looked at this know from this kind of academic perspective and it might seem like we went into the nitty-gritty details I want to kind of leave you with one one small idea of how you can maybe extend this concept right if you think about what we did tonight is we linearize the system about one operating point right about trim steady state straight and level cruise at eighty five meters per second right now you might have a system which can operate at a lot of different speeds right and a lot of different altitudes a lot of different attitudes a lot of different orientations basically a lot of different trim points right so the linear model that we obtain tonight is sort of only good at one of these points like so for example you might see an operational envelope of an aircraft you know spread out in two dimensions like a long mock and altitude right and what we did tonight is we picked you know some point in this operating envelope leg you know maybe here here's where what we did tonight right here is eighty five meters a second at some altitude or whatnot right and we got ourselves an A and a B matrix about this point so now that we have this linear system we could go ahead and do control design we could do stability analysis we could do whatever we wanted about this location in fact that's what we're gonna do in one of our future videos is the immediately following video is we're gonna start talking about how do I analyze the system about this operating point what can I learn about the aircraft should behavior that system power at that trim point right well if you move the trim point of the vehicle like you start flying at some other location the system is it might behave differently you might have to redo this analysis that we just talked about and generate another set maybe let's call this like an A two and a B to a different set of linear models at that location and redo all of your analysis and you can imagine doing this at different locations in the envelope in other words you would have different linear models for different points in your operating regime right so this opens up a large variety of additional analysis of what we can do in the system but now we've built the tools to enable us to go ahead and linearize the system about any of these points in in in state space right where you'd like to look at it so I think that's pretty exciting and I hope you do too so that being said I hope you enjoyed the discussion tonight and if so I also hope you'll consider subscribing to the channel surprisingly again if you just scroll little ways down and click on that subscribe button it really does help me continue making these videos and if you liked it or if you didn't like the video tonight please leave a comment below and let me know I'd love to get some feedback and also we can get a discussion going among like-minded individuals who are also studying this topic so until I see you at one of these future videos I think I'll sign off for now so I'll talk to you later bye

Channel: Christopher Lum

Views: 6,394

Keywords: Linearization, linearize, numerical linearization, linearize perturbation method

Id: 1VmeijdM1qs

Length: 104min 45sec (6285 seconds)

Published: Thu May 21 2020