tag:blogger.com,1999:blog-6476989946886057900.post6739532513963428656..comments2019-12-19T15:25:37.909-08:00Comments on <small> <i> the </i> </small> Mechanical Money <small> <i> blog </i> </small>: For Tom: The GDP ObjectRoger Sparkshttp://www.blogger.com/profile/01734503500078064208noreply@blogger.comBlogger38125tag:blogger.com,1999:blog-6476989946886057900.post-87553415959225165252016-04-05T14:59:29.553-07:002016-04-05T14:59:29.553-07:00"I wonder if Excel has a iteration loop for a..."I wonder if Excel has a iteration loop for a single cell, that could run and continually sample a "future" cell?"<br /><br />You can turn something like that on here:<br /><br />Excel->options->formulas->'Enable iterative calculation'<br /><br />There's also what they call their "Solver" (really a single root finder) and the fancier "non-linear equation solver" which you need to install their analysis package for. The analysis package has a fancier linear regression capability (than the usual one in their formulas, or built into their plots), and it has matrix manipulation capabilities like matrix multiplies and matrix inverse. Octave is much better for matrix manipulations though!Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-23262998841525775822016-04-05T14:07:29.281-07:002016-04-05T14:07:29.281-07:00Tom, Just thinking here. I'm reading your mate...Tom, Just thinking here. I'm reading your material (and more) about control systems. I can now see the connection between SIM and control loops. So yes, I think Eq. 10 (the Master Equation) would be a control object (programming concept of an object). <br /><br />Any feedback into the control object would necessarily be from the past (in a time sense) so it could be added to any of the four controlling parameters. [I wonder if Excel has a iteration loop for a single cell, that could run and continually sample a "future" cell?] <br /><br />I think that from a system control viewpoint, the Master Equation is an integrator summing four different inputs including one sample from memory. The G input is master. and inputs 2 and 3 are successive and diminishing in effectiveness.<br /><br />More inputs could be added at any of the entry points.<br /><br />I am starting to get a glimmer of how the Z transform works (I think).<br /><br />This is all very cool. Thanks! Roger Sparkshttps://www.blogger.com/profile/01734503500078064208noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-21683454372632955592016-04-05T10:16:20.296-07:002016-04-05T10:16:20.296-07:00It's unfortunate that the star (*) is the norm...It's unfortunate that the star (*) is the normal symbol for convolution, since it's used to represent multiplication in text equations and computer programs. But Wikipedia has a good article on it (complete with animations):<br /><br />https://en.wikipedia.org/wiki/Convolution<br /><br />It's easier to visualize in discrete time: basically everywhere an input signal to a system is non-zero, then for the output, replace that with a scaled and shifted system impulse response, and then add all those scaled and shifted impulse responses together: that's the superposition concept at play: you can't do that with a non-linear system where superposition doesn't hold. (it even works for a time-varying linear system, only the impulse response in that case is a function of time: it's not fixed like for an LTI system).<br /><br />Tables of Laplace and Z-transforms (and their inverses) make solving convolution integrals (convolution sums in discrete time) unnecessary if you can describe the impulse response and the system input as shifted and weighted sums of common functions in the table.<br /><br />A classic case is the action of a digital to analog converter (DAC). A DAC transitions you from the discrete time world to the continuous time world. Typically they work by a "sample and hold" or zero order hold (ZOH) technique, meaning they hold their output voltage constant during each sample period. You can think of the DAC's impulse response then as a shifted rectangle (which in turn are two step functions, one up and one down, added together: u(t) - u(t-T), where u(t) here is the unit step function). Essentially that impulse response is being convolved with the discrete time input to the DAC. In continuous time a discrete time signal's spectrum is repeated every factor of Fs (the sampling frequency). But the frequency response of a rectangular pulse is a modified sinc() function (having the general form sin(a*f)/(a*f)). So in the frequency domain, you multiply the spectrums (since you're convolving in the time domain), thus the spectral images of the discrete time signal are attenuated as you go further out in frequency in the continuous time output of the DAC. Also a shift in one domain is a phase change in the other. So a time shifted rectangle has a frequency response which is a phase shifted version of the spectrum of a 0 centered rectangle: the magnitude doesn't change. Another principal is that signals which are wide in the time domain, are narrow in the frequency domain, and vice versa. So a wide rectangle (a long sample period) means a narrow sinc in the frequency domain. A Gaussian bell curve is the same in both domains other than its width. So a wide Gaussian in one domain transforms to a narrow Gaussian in the other.<br /><br />Still, a DAC will produce many undesirable (though attenuated) harmonics of the fundamental (the desired signal), and you generally add an analog filter to its output to eliminate the ones you don't want.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-31700083861850711422016-04-05T10:07:58.933-07:002016-04-05T10:07:58.933-07:00Tom, I just posted A Master Equation for SIM Model...Tom, I just posted <a href="http://mechanicalmoney.blogspot.com/2016/04/a-master-equation-for-sim-models-using.html" rel="nofollow">A Master Equation for SIM Models Using the GDP Object</a> . I labeled it draft and welcome comments for improvement. I think I am satisfied with this but are you? My understanding kept evolving as I wrote.<br /><br />I will look at your comments next.Roger Sparkshttps://www.blogger.com/profile/01734503500078064208noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-84583835648216500142016-04-05T08:23:48.851-07:002016-04-05T08:23:48.851-07:00Above, where I wrote:
"1st order linear MATR...Above, where I wrote:<br /><br />"1st order linear MATRIX differential equation"<br /><br />I should have written:<br /><br />"1st order linear VECTOR differential equation"<br /><br />The vector being x. It needs a square matrix A to fill out the equation, but A is not the state. A matrix differential equation would be a case where a matrix X was involved:<br /><br />X' = A*X<br /><br />Which is really just a bunch of vector differential equations:<br /><br />x1' = A*x1<br />x2' = A*x2<br /><br />etc, where xi is the ith column of X. You can write the optimal, time varying, continuous time Kalman filter gain equation in a matrix differential equation form, but unfortunately it's non-linear... I forget what it is, but it's something like:<br /><br />X' = R*X + X*R + X*Q*X + F<br /><br />or something horrible like that, where X is the Kalman gain. (The Kalman gain is a time varying gain you need to perform an optimal state estimate). It's the Riccati equation I mentioned once, only in a matrix differential form.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-89236492413274938932016-04-05T07:47:36.154-07:002016-04-05T07:47:36.154-07:00In case you're wondering, the analog of an int...In case you're wondering, the analog of an integrator in the discrete world is an accumulator. So instead of this:<br /><br />x' = u<br /><br />You have:<br /><br />x[n+1] = x[n] + u[n]<br /><br />Taking the Laplace transform of the 1st:<br /><br />X(s)*s = U(s)<br /><br />So that's why an integrator is 1/s: X/U = 1/s<br /><br />In the discrete time world, what do we have? Well the way I wrote it, x[n+1] is "delayed" by -1 (i.e. advanced by one) so:<br /><br />X(z)*z = X(z) + U(z)<br /><br />X/U = 1/(z-1) = z^-1/(1-z^-1)Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-44063772935951165722016-04-05T07:39:33.062-07:002016-04-05T07:39:33.062-07:00One last thought: using the state space representa...One last thought: using the state space representation allows you to write the system differential equation as a 1st order linear MATRIX differential equation:<br /><br />x' = A*x<br /><br />Just 1st order!... but in terms of a scalar differential equation, it's not 1st order, it's order is the dimension of x! So if x is Nx1, the it'll be Nth order. All those terms in the characteristic equation (a*s^n + b*s^(n-1) ... etc, down to c) are all just differentiated x's... in a differential equation they'd be terms like a*x''' + b*x'' + ... + c = 0. So taking the Laplace transform of a differential equations replaced nth derivatives with s^n, and it replaces integrations over time with factors of 1/s, double integrals with 1/s^2, etc.<br /><br />All the same applies for Z transforms, except that instead of differentials and integrals you look for cases where the time index is different. If a term is delayed by one x[n-1] then in the Z domain it's X(z)*z^-1. If it's delayed by 2 then it's X(z)*z^-2, etc. A delay in digital circuits is often marked with Z^-n, where n is the magnitude of the delay.<br /><br />Laplace transforms and Z-transforms can make life easy once you get the hang of it. You can look up in a table what the transforms are for many common functions, and if you know a few rules, you can adapt those to your circumstance.<br /><br />In the end with Z transforms, it boils down to the same thing for SISO systems: a ratio of polynomials in z. So it has a root locus, just like its continuous time cousin in s.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-19511172038380261602016-04-05T07:22:08.383-07:002016-04-05T07:22:08.383-07:00Then he factors P(s) into it's numerator and d...Then he factors P(s) into it's numerator and denominator... so you can see what K does to the roots of the closed loop system. I forget what he calls them, but lets call them N(s) and D(s). Then K the gain moves the roots of the closed loop system from those of D(s) (when K = 0) to those of N(s) (when K = infinity). The zeros of a system (open loop, those are the roots of N(s)) have an effect on a system, but they don't determine it's stability or it's "speed" or oscillatory behavior... those critical things are all determined by the roots of D(s) (again, open loop, not closed loop), which again are the roots of<br /> <br />det(A - s*I)<br /><br />which are the eigenvalues of A, which are the poles of the system, which are the roots of characteristic equation... a lot of names for the same thing!<br /><br />But the interesting part is how as K goes from 0 to infinity the poles of the *CLOSED LOOP* system move from D(s) to N(s). The poles move to the zeros.<br /><br />Now imagine that K isn't just a constant gain, but also has a form like this:<br /><br />K(s) = Nk(s)/Dk(s)<br /><br />It's easy to work out (in terms of numerator polynomials and denominator polynomials) what the closed loop system is. It's best to just reduce it to one numerator polynomial divided by one denominator polynomial so it's easy to see (IMO).<br /><br />As he shows, K a constant doesn't always do the trick in terms of being able to stabilize a naturally unstable system. Sometimes you have to add dynamics to K as well. In the one case I saw, he added a zero (which is called a differentiator), so K had a form like<br /><br /><br />K(s) = a + b*s<br /><br />Or what's known as proportional (a) differential (b) control: a PD controller. If you add in an integrator as well (PID), it's something like a + b*s + c/s.<br /><br />1/s is the Laplace transform of an integrator, and s is the Laplace transform of a differentiator.<br /><br />The problem with differentiators is they tend to amplify noise (they are high pass filters), whereas integrators tend to reduce noise (low pass filters)... but noise may be less important than stability!<br /><br />Sometimes what happens with your controller (K) is that it actually has a model of the plant built right into it... so K has within it P(s)... now why would that be? Well that would be to estimate ALL the states of P... recall, we can express P's polynomial as a set of states with x' = A*x.<br /><br />That's precisely what a Kalman filter is: it's a little model of the system P which estimates all the states. It's married together with a control law which then tries to modify the states (change the characteristic equation). If you can accurately model all the states, then you can put the poles of the system wherever you'd like.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-28879031744434114132016-04-05T07:21:17.503-07:002016-04-05T07:21:17.503-07:00Roger, you're not slow to learn!!! Shoot, I...Roger, you're not slow to learn!!! Shoot, I've been drowning you in stuff it took me a lot of time to learn.<br /><br />Now Brian didn't show in that 1st video (perhaps earlier ones he does) how gets the equations for those "closed loop" transfer functions in the 1st place. It's actually very simple. You can do the exact same thing with Z transforms too BTW. So open loop you may have:<br /><br />Y(s) = P(s)*U(s)<br /><br />Where y(t) is the output and Y(s) is it's Laplace transform, and same for u(t) and U(s). P(s) is the transform for the "plant" (the system, like SIM is). In the time domain, the "impulse response" of the plant in time (call it p(t)... like what happens if you hit SIM with just $1 of gov spending in the first period and then back to $0 afterwards) you'd say:<br /><br />y(t) = p(t)*u(t)<br /><br />too, but the "*" symbol means "convolve with" rather than multiply in that case. "convolve with" means for every input in u(t) you put a copy of p(t)'s impulse response. The impulse response times 20 is what you'd get if you didn't accumulate at the end of Jason's version of SIM, because he just has G at $20 for period 1, and $0 after.<br /><br />So multiplication in the frequency domain is convolution in the time domain. The system is linear, so to find the response to ANY input is just the sum of impulses responses to a bunch of impulses added together.<br /><br />Now the converse is true as well: convolution in the frequency domain, is equal to multiplication in the time domain... so if two signals are multiplied in the time domain (like what a mixer does in radios), then in the frequency domain their spectrums are convolved.<br /><br />But back to Brians videos:<br /><br />Y(s) = P(s)*U(s), so the transfer function is just (dropping the s argument):<br /><br />Y/U = P<br /><br />And because most realizable systems have Laplace transforms that are ratios of polynomials, those are just all polynomials in s. Now, these are what are called SISO systems (single input single output). Translating this into a "state space" representation still means that A can be NxN in general (lots of states = lots of terms in the polynomials), but it means that the y output and the u input are both just scalars. Otherwise you have LOTS of impulse responses: from each u to each y. It's likely however, that they will all have the same polynomial in the denominator: only the numerator polynomial will change (most of the time). That denominator polynomial (in s) is found from the characteristic equation of A:<br /><br />det(A - s*I) = 0<br /><br />But Brian starts with just the polynomials themselves, which is fine. <br /><br />The other thing you need to know is that addition in the time domain is addition in the frequency domain. So when you take Y and feed it back to the input by subtracting it from U (negative feedback... the usual) you get after "closing the loop":<br /><br />Y = P*(U - Y)<br /><br />or<br /><br />Y + P*Y = P*U<br />Y(1+P) = P*U<br />Y = P*U/(1+P)<br />Y/U = P/(1+P)<br /><br />Now he puts a gain K in there after in the forward path after the summing junction, so it's <br /><br />Y/U = K*P/(1-K*P)<br /><br />Where that "1" is really a 1 and not an I (not an identity matrix in this case), because of course it's a SISO system. These functions of s are all just scalars.<br /><br />Now K can also be a function of s, so it's K(s), but the easiest thing to do is just make it a scalar K. And convolution by a constant K is the same as multiplication by a constant K.<br /><br />Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-70028340003979284702016-04-05T05:49:23.318-07:002016-04-05T05:49:23.318-07:00Tom, Thanks for sharing this link to Brian Douglas...Tom, Thanks for sharing this link to Brian Douglas' video series. Very interesting. <br /><br />Sorry that I am so slow to learn. It seems like I need to 'internalize' things, making them part of my mechanical system, the way I perceive things.<br /><br />That said, I was thinking about matrix's last night. What a beautiful system of organization. Of course, all they are is a regular system of identification. Each point on the grid corresponds to a unique designation. You understand, I am just opening my eyes to the possibilities! Thanks again for pointing the way.<br /><br />I am slowly beginning to understanding the link between SIM and control systems. I recognized one of the equations (a part) that Brian Douglas displayed, so you will (hopefully) like the new approach to my SIM equations. Hopefully, a post today, and more SIM Excel examples later in the week.Roger Sparkshttps://www.blogger.com/profile/01734503500078064208noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-60377185557628848292016-04-04T21:53:13.658-07:002016-04-04T21:53:13.658-07:00With the array in it? You mean [A B; C D]? Don'...With the array in it? You mean [A B; C D]? Don't worry about that... it doesn't do anything. You can comment it out. It shouldn't cause anything to hang though.<br /><br />Check it out Roger: you probably won't follow all this (I didn't either actually... I used to draw these things myself, but it's been a while). Anyway, this guy has two videos that are kind of fun to watch: he shows you what happens to the roots of a system as you change the control function. Here's the 1st video:<br /><a href="https://www.youtube.com/watch?v=eTVddYCeiKI&feature=iv&src_vid=jb_FiP5tKig&annotation_id=annotation_831977" rel="nofollow">https://www.youtube.com/watch?v=eTVddYCeiKI&feature=iv&src_vid=jb_FiP5tKig&annotation_id=annotation_831977</a><br /><br />It's all actually very geometrical. His "rules" may seem mysterious, but it all comes down the the geometry of the poles and zeros in the plane: their distance and their angles. Anyway, it's kind of fun, because he's got the intuition nailed, so he can draw them fast.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-21445332687429046502016-04-04T21:43:47.437-07:002016-04-04T21:43:47.437-07:00Tom, I'm deep into a rewrite. I think this eff...Tom, I'm deep into a rewrite. I think this effort is better. Hopefully, I will post in draft form tomorrow. <br /><br />I tried Octave with the array in. It hangs someplace after figure one. I will play with it more. Next, I will read the instructions. So far, I am just going on intuition based on some C++ programming I have done, and on the comments and wonderful code that you wrote.Roger Sparkshttps://www.blogger.com/profile/01734503500078064208noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-25418952770512610332016-04-04T20:20:31.244-07:002016-04-04T20:20:31.244-07:00Roger, you're right: there's no "box&...Roger, you're right: there's no "box" marker: it's a square. From Matlab help:<br /><br /> b blue . point - solid<br /> g green o circle : dotted<br /> r red x x-mark -. dashdot <br /> c cyan + plus -- dashed <br /> m magenta * star (none) no line<br /> y yellow s square<br /> k black d diamond<br /> w white v triangle (down)<br /> ^ triangle (up)<br /> < triangle (left)<br /> > triangle (right)<br /> p pentagram<br /> h hexagramTom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-78593902958872393552016-04-04T13:40:20.178-07:002016-04-04T13:40:20.178-07:00... another fun fact, even though 3rd order and hi...... another fun fact, even though 3rd order and higher solutions are notoriously difficult to solve by hand, with a few geometric rules in the S or Z plane, you can still approximately sketch out their "root locus" diagrams as a function of the feedback gain in a system. It's pretty cool...Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-77662532979363896872016-04-04T13:37:14.629-07:002016-04-04T13:37:14.629-07:00specifically, the "characteristic equation&qu...specifically, the "characteristic equation" is the polynomial equation resulting from this:<br /><br />det(A - lambda*I) = 0<br /><br />The values of lambda that satisfy that are the roots of the characteristic equation and the poles of your system.<br /><br />That reduces to a regular old scalar coefficient polynomial... so if A ia a 2x2, it's a quadratic formula. "det" is the determinant. Even calculating the determinant gets to be a major pain if you do it symbolically with a 3x3 or a 4x4, etc. Except in special cases. There's never a reason to do it by hand really... except on exams!<br /><br />Fun fact: If you replace the scalar lambda in the characteristic polynomial, with the *matrix* A, then A satisfies it! (provided, of course, you replace the 0 with a matrix of zeros, and interpret A^n as (n-1) matrix multiplies of A with itself). In other words, A satisfies it's own characteristic equation. Strange, eh?Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-56050823958544232182016-04-04T13:26:53.274-07:002016-04-04T13:26:53.274-07:00If you put a semicolon at the end, the line won...If you put a semicolon at the end, the line won't print. So <br />[A B; C D]<br /><br /><br />didn't actually do anything except form that one big matrix, but it wasn't assigned to anything, so it has no function except to print to your console. The thing is, it probably makes a mess, so I can understand you not wanting to see it! :D<br /><br />for a cleaner way to print to console look up "disp" and "display" ... I was just too lazy to use those.<br /><br />Matrix math... one of the Z and S plane links comes from a page all about mechanical systems: it turns out the RLC circuit Jason was using for an analogy is also analogous to a simple 2nd order mechanical analogy: a mass (usually pictures as a horizontal cart with frictionless wheels & bearings), a massless spring, and a "dashpot" (like a dampening device). The have the diagram right there on one of those pages. It fits in well with your "mechanical perspective" I think.<br /><br />My "feedback control" script above doesn't use any matrices (just scalars) so it's easier to follow. It should be clear in that case how adding feedback move the system's pole, but you can run it and verify that for yourself. Since SIM is a one pole system, it's directly analogous, especially when people start moving the pole around with gamma or alpha2 (as Jason and Ramanan do, resp.).<br /><br />BTW, here's that page with the mass, spring and dashpot... and a servo motor actuator as well (to provide a forcing function u, much like Jason's voltage generator in his RLC circuit):<br />http://nupet.daelt.ct.utfpr.edu.br/_ontomos/paginas/AMESim4.2.0/demo/Misc/la/SecondOrder/SecondOrder.htm<br /><br />2nd order systems are very popular on exams because they're just complex enough to be interesting, yet you can solve them with nothing more complex than the quadratic formula (to find the roots of the system). Solving 3rd order polynomial problems and higher was not generally a fun or fruitful project with your pocket calculator back when I was in school! And of course, in general there's no closed form solution for 5th order or above (except special cases) as Galois proved more than a century and a half ago. (If you've ever seen the closed form for a 4th order polynomial, it's UGLY!). Matlab easily solves all that numerically though... that's essentially what eig() does.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-18594700778311719882016-04-04T08:15:49.687-07:002016-04-04T08:15:49.687-07:00Tom, Those Z and S plane figures (that you referen...Tom, Those Z and S plane figures (that you referenced) are very interesting. They invite careful examination and a lot of study.<br /><br />OH, I changed the plot style/color from rb to sr. Matlab "red box' to Octave "square red". Now I get 10 scattered squares. <br /><br />Also changed [AB; CD] to [AB; CD]; . Now I also have two curve displays. The curves are not stable and equal replication-to-replication so maybe more needs to be done.<br /><br />It astonishes me that you wrote the code nearly error free the very first time. Shows great skill!<br /><br />I am trying a new approach on "The GDP Object". My argument must convince folks skilled in math, as well as convincing me. That is a high standard, not yet met.<br /><br />OH, more. I spend some time working farther into matrix math. Yes, I see how it works. Now I need to get a better skill level in using the tool.Roger Sparkshttps://www.blogger.com/profile/01734503500078064208noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-20533773556089503502016-04-03T19:35:53.329-07:002016-04-03T19:35:53.329-07:00In the script, this:
Y = (1,N);
Was supposed to ...In the script, this:<br /><br />Y = (1,N);<br /><br />Was supposed to be this:<br /><br />Y = zeros(1,N);<br /><br />Also instead of<br /><br />L = 0.51;<br /><br />It should be<br /><br />L = -0.51;Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-86367328207329857622016-04-03T19:28:54.320-07:002016-04-03T19:28:54.320-07:00OK, a simple feedback control example: imagine the...OK, a simple feedback control example: imagine the following single state system:<br /><br />x[n+1] = a*x[n] + b*u[n] + g*w[n]<br />y[n] = c*x[n] + m*v[n]<br /><br />Where c = 1, g = 0.3, m = 0.3, b = 1 and a = 1.01, x0 = 0 and u[n] = 1, n >= 0, else 0, and w and v are N(0,1) (zero mean, unity variance additive white Gaussian noise (AWGN).<br /><br />Then xhat = y is an unbiased estimate of x. The pole of this system is simply a, which is just outside the unit circle, so it should be unstable. But now lets apply a gain L to xhat and feed it back:<br /><br />x[n+1] = a*x[n] + b*u[n] + L*x[n] <br /><br />Where I've eliminated the noise terms for simplicity. So as long as L is on (-0.01,-1.01) we can stabilize the system:<br /><br />X(z) = (a+L)*X(z)*z^-1 + b*U(z)<br /><br />X(z)/U(z) = 1/(1-(a+L)*z^-1) = z/(z - (a+L)), so the pole is at (a+L). If we don't want it to oscillate, and we want to dampen the control signal u, then a good choice would be a+L on (0,1). Let's put the pole at 0.5, which means L = 0.51.<br /><br />a = 1.01;<br />b = 1;<br />L = 0.51;<br />g = 0.1; % "plant noise" standard deviation<br />m = 0.1; % measurement noise standard deviation<br />u = 1;<br />x0 = 0;<br />N = 100; % samples<br /><br />Y = (1,N);<br /><br />% first with no feedback control compensation:<br />x = x0;<br />for k=1:N<br /> Y(k) = x + m*randn<br /> x = a*x + b*u + g*randn<br />end<br /><br />figure(1),plot(Y),title('No feedback control');<br /><br />% Now with feedback control compensation:<br />x = x0;<br />for k=1:N<br /> Y(k) = x + m*randn<br /> x = a*x + b*u + g*randn + L*Y(k);<br />end<br /><br />figure(2),plot(Y),title('With feedback control');<br /><br />% END script<br /><br />With more complex systems, you can determine what the poles will do ahead of time with different gains L when you "close the loop." Generally you make a plot of how they move as function of the gain. And sometimes your "gain" has states of it's own. <br /><br /><a href="http://12000.org/my_courses/UCI_COURSES/CREDIT_COURSES/spring_2005/MAE_171_UCI_spring_2005/HWs/HW7/HTML/graphics/hw7__42.gif" rel="nofollow">Here's a typical diagram</a> in the Z-plane showing how two poles move as a function of the feedback gain in relation to the unit circle (black). You can see (red lines) how they (the poles) move along the real line until they meet, and then they branch off into conjugate complex numbers, circle around, and meet again on the other side of zero, where they again split from each other on the negative real line. Any one of those (pair of) positions on the red curves corresponds to a particular gain. So the designer would choose a gain that puts them where he's like.<br /><br /><a href="http://nupet.daelt.ct.utfpr.edu.br/_ontomos/paginas/AMESim4.2.0/demo/Misc/la/SecondOrder/SecondOrder_fichiers/image070.png" rel="nofollow">I thought this was a nice diagram too</a> showing what the impulse response of poles look like at different locations around the S-plane (so continuous time this time). Now instead of the unit disk (the unit circle and it's insides) the stable region is the right half plane (and the imaginary axis). <a href="http://dspcan.homestead.com/files/Ztran/s-z-mapa.gif" rel="nofollow">This diagram</a> provides a nice map between the S-plane and the Z-plane. <a href="http://dspcan.homestead.com/files/Ztran/zlap-z.gif" rel="nofollow">This one's good too.</a>Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-58297077925890690852016-04-03T18:26:23.804-07:002016-04-03T18:26:23.804-07:00Also, Bill makes an interesting comment here.Also, <a href="http://informationtransfereconomics.blogspot.com/2016/03/an-rlc-circuit-with-r-s-and-l-f.html?showComment=1459707901983#c5380084260083725576" rel="nofollow">Bill makes an interesting comment here</a>.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-52329863314108041512016-04-03T14:28:18.805-07:002016-04-03T14:28:18.805-07:00I just noticed an error with my expression for y a...I just noticed an error with my expression for y above (not in the script). I wrote:<br /><br />y = C*x + B*U(:,k) + M*randn(Ny,1);<br /><br />But it should be<br /><br />y = C*x + D*U(:,k) + M*randn(Ny,1);<br /><br />and of course in term of the script, where we're recording every y, then it's<br /><br />Y(:,k) = C*x + D*U(:,k) + M*randn(Ny,1);<br /><br />It turns out SIM is a bit of an odd beast in this regard... it actually has a non-zero D. If you think about it for a system like a rocket, servo motor or inverted pendulum that you're trying to control, you don't normally have a non-zero D: that would just be taking a linear combination of the input control signal u (which you presumably already know) and adding that to your measurements. Why would you do that?<br /><br />And like I say above, you're lucky if you know A, even approximately, but as a control engineer you can't usually have any say over what A is: it describes the system you've been instructed to control. You might be able to determine that the system is uncontrollable (depends on A and B) or unobservable (depends on A and C)... which could lead to a redesign of the system, or tell something about the stability of states for which there's no controllability or observability (you'd better hope in both cases they are at least stable). In the same manner, your analytical feedback to your customer could lead to adopting different or more sensors (C) or adding more controls/actuators (B), or improving the quality of components to decrease or alter the character of the (thermal) noise (M or G), ... or giving up on the project altogether!Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-4328781553413143522016-04-03T13:53:03.073-07:002016-04-03T13:53:03.073-07:00Roger, if you want to see an application of this l...Roger, if you want to see an application of this linear system theory stuff, and particularly a "root locus" (i.e. predicting system behavior based on the location of the systems roots (poles)), Jason tried setting gamma = 10 and noticed that the system oscillated, and thought that was interesting. In terms of how gamma affects the location of the system's pole (just the scalar A in the case of SIM), and how that in turn can explain and predict SIM's behavior, I left a <a href="http://informationtransfereconomics.blogspot.com/2016/03/an-rlc-circuit-with-r-s-and-l-f.html?showComment=1459712467572#c1832579359314853328" rel="nofollow">couple more comments in response on his RLC circuit post</a>.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-13619668876865660842016-04-03T11:25:31.597-07:002016-04-03T11:25:31.597-07:00... economics is a classic case requiring "ro...... economics is a classic case requiring "robust control": even if the system were linearizable (which it may not be), you don't really know A, B, C, D, G or M... so you have to build a feedback control law (and state estimator) that works (to some minimal level of performance: stability at least for example) for a whole set of possible plants! It's not an easy chore. Perhaps impossible. You've got to hand it to Taylor for giving it a stab anyway.<br /><br />Or you can look at it as a non-linear control problem... even harder!<br /><br />In fact, just figuring out a nominal A and the other matrices and the range of uncertainty on them (system identification) is proving to be pretty much impossible I think.<br /><br />BTW, you can handle "colored" noise (not white: white meaning the autocorrelation function is an impulse... i.e. all frequencies represented equally) with this setup by adding extra states designed to color the noise.<br /><br />Also, some economists think the system is stable to start out, and what makes it unstable are attempts to add feedback. So there's that too. Others think it's not stable, but you can't improve it any with feedback, so it's no use trying.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-78737948048323176642016-04-03T11:11:40.237-07:002016-04-03T11:11:40.237-07:00... yeah, I imagine [A B; C D] would make a mess w...... yeah, I imagine [A B; C D] would make a mess with 10 states!... I originally had 3, but I bumped it up so you'd be pretty likely to have some unstable poles outside the unit circle. I have no idea what your output looks like: probably just a mess with a random input. The impulse response (or a step response) should look pretty clean, but you probably get gigantic numbers! 100 iterations on one or more unstable poles can send it to infinity fast. You might have to cut down on the number of sample times to avoid overflows, even with double precision numbers.<br /><br />Now we could add noise to the plant with covariance G'*G, and noise to the measurements with covariance M'*M:<br />G = randn(Nx,Nx);<br />M = randn(Ny,Ny);<br />% the "random" u is already like a noise with<br />% covariance eye(Nx) (i.e. an Nx sized identity matrix)<br />...<br /><br />y = C*x + B*U(:,k) + M*randn(Ny,1);<br />x = A*x + B*U(:,k) + G*randn(Nx,1);<br /><br /><br />and then we could build a Kalman filter to optimally track all 10 states and plot the principal axis length of the error covariance based on just those two outputs... but we can leave that for another day (like maybe never, since you've probably had enough!... Lol).<br /><br />Fun fact about the function eig(): it knows how to diagonalize a square matrix, but that's it. Of course not all matrices are square, but even the square ones aren't all diagonalizable. The diagonal values being the eigenvalues of course. Once simple one that isn't (which comes up all the time) is the A matrix associated with the states position, velocity, acceleration, rate of acceleration (jerk), etc. Such a simple matrix and system, but one of the extremely rare cases that's not diagonalizable. With just position and velocity, in continuous time that's:<br /><br />A = [0 1; 0 0]; % Not invertible or diagonalizable!<br /><br />And in discrete time (invertible, but not diagonalizable):<br /><br />A = [1 T; 0 T]; % where T is the time step.<br /><br />If you try eig(A) in either case I'd guess it'll fail. That doesn't mean there are no eigenvalues... the problem is there's not enough eigenvectors. But here's the thing: in either case, A is just a tiny epsilon away from a diagonalizable matrix:<br /><br />A = [1 T; 1e-10 1]<br /><br />Now you can diagonalize A just fine. There's no such problem with singular values... but then they don't do EVERYTHING, so you need the eigenvalues in many cases.<br /><br />Once you build a state estimator, you'll have a state estimate xhat, and then you can feed that back to stabilize the system instead of manually moving the poles like I did:<br /><br />x = A*x + B*U(:,k) + G*randn(Nx,1) + L*xhat;<br /><br />Where L is your "control law"... the idea being in the real world you don't have access to A (other than a model of it): you only have y and you know what u is hopefully. So you form xhat from y and u, and use it to stabilize the system.<br /><br />That's how they can stabilize that inverted pendulum. They know the angle of each segment (with sensors)... and they can do a transformation on that to move the cart to always restore it to "equilibrium." The feedback moves all the poles inside the unit circle. I'm only talking about when it's already close to equilibrium... I have no idea how they move the cart to stack the segments up.<br /><br />A potential economics application? You have an unstable system (or one with undesireable dynamics anyway) you'd like to stabilize/make-more-desireable perhaps? I'm sure I'm only the billionth person to have thought that! Lol. (That's what the Taylor rule is all about, BTW: that's L the control law).<br /><br />Here's another application: a servo motor. You might be able to stabilize the system so that a given voltage corresponds to a certain angle of the drive shaft... but stable is not always good enough. Perhaps it oscillates a bit when it gets there... the oscillations die down to zero, but you may not be able to tolerate any. So you can use feedback to dampen those out, or get rid of them completely.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.comtag:blogger.com,1999:blog-6476989946886057900.post-12584067854692362482016-04-03T11:02:10.410-07:002016-04-03T11:02:10.410-07:00This comment has been removed by the author.Tom Brownhttps://www.blogger.com/profile/17654184190478330946noreply@blogger.com