Sunday, April 10, 2016

SIM9: highlighting the difference between canonical and non-canonical and between average and instantaneous rates

Update: simple circuit circuit analogy (and interactive simulation) here.

The math for these plots will be posted at a later time. Please refer to SIM8 for now. A couple of quick notes: the canonical output y from SIM8 has been replaced here with w1. The curve w2 is new. Both w1 and w2 are time shifted with respect to their corresponding SIM outputs. I've only included the tax collection rate T (so w1 = T) in the plots for an output (which I've added a subscript of "1" to), but I've worked out the math for the full set of outputs (Y, T, YD and C), which are all time shifted (along with the input G) with respect to the canonical form. The state variable H is aligned in time with the canonical state variable x. But for the others we have:

G[n] = u[n-1]
T1[n] = w1[n-1]

w1 is a time shifted version of the average tax collection rate over the coming period. w2 is the time shifted instantaneous tax collection rate at the start of the period. This instantaneous rate I call T2, and thus following the above we have T2[n] = w2[n-1]. All the outputs will have two versions like this: average and instantaneous. The older average versions will use the same C and D system matrices, but the newer instantaneous ones will need new C and D matrices. Those details are coming in later updates to this post.

I've reluctantly decided to abandon y for the name of the canonical output variable in favor of w. I'm tired of the name collisions between SIM variables and canonical variable names, so I'll try to sort this out. Also, the lower case Latin letter versions with the "dots" are OK for the average rates, but I found it easier to just adopt a whole new set for continuous time based on lower case Greek letters. I'll explain the convention in later updates.

This difference between the average and the instantaneous and between the canonical and non-canonical variables has been bothering me for some time, so I think will finally completely sort all this out, but I need wait until further updates to write it down. I thought I'd put the plots up now however: perhaps they can provide a little bit of clarity (probably not! Lol):

UPDATE: 2016.04.12: Updated the figures and added a second canonical one with the same curves but different annotations.

Figure 1: Canonical #1
Figure 2: Canonical #2 (curves are the same as Canonical #1, but annotation is different)
Figure 3: Non-canonical
Figure 4: All curves
Update 2016.04.18: if we define P (previously referred to as $X$ in SIM7) as:
$$P = 1-\alpha_1 (1-\theta) \tag 1$$
Then we can write the elements of the "instantaneous" version of the state space representation measurement matrices ($C_{inst}$ and $D_{inst}$ in terms of their 1-period-sliding-window-average state space measurement matrices ($C_m$ (or $C_{avg}$) and $D_m$ (or $D_{avg}$)) as follows (where $i$ is the index of the $i^{th}$ row (i.e. $i^{th}$ element) of these 4x1 matrices):
$$C_{{inst}_i} = -\frac{C_{{avg}_i} P}{\alpha_2 \theta} a \tag 2$$
$$D_{{inst}_i} = D_{{avg}_i} - \frac{C_{{avg}_i} P}{\alpha_2 \theta} (b - B) \tag 3$$
Up to this time we've only calculated the $C_m$ and $D_m$ matrices. Here are the full set of values (including for the op-amp implementation of SIM) for the default parameter values of SIM (i.e. $\alpha_1 = 0.6,\, alpha_2 = 0.4$ and $\theta = 0.2$.

Name value 1/value Instantaneous  1/Instantaneous
A 0.846153846 N/A (discrete time)
B 0.615384615 N/A (discrete time)
a_cont -0.167054085 -5.9860852970
b_cont 0.668216339 1.4965213242
Cavg = Cm Cinst 
CY 0.769230769 1.3000000000 0.835270423 1.197217059
CT 0.153846154 6.5000000000 0.167054085 5.986085297
CYD 0.615384615 1.6250000000 0.668216339 1.496521324
CC 0.769230769 1.3000000000 0.835270423 1.197217059
Davg = Dm Dinst 
DY 1.923076923 0.5200000000 1.658918307 0.602802438
DT 0.384615385 2.6000000000 0.331783661 3.014012191
DYD 1.538461538 0.6500000000 1.327134645 0.753503048
DC 0.923076923 1.0833333333 0.658918307 1.517638818

A preliminary implementation of this instantaneous version of the op amp SIM circuit is here.

Same functionality with only 5 instead of 7 op amps.

Same functionality, but preserving H volts on capacitor and using 6 instead of 7 op amps.

Same thing, but with integrator on the Y output (1/2 of a ping pong integrator)

Same thing, but with full ping-pong integrator and synchronized start switch.  

Note to self: to-do list (2016.04.10: 1:17 PM).
Updated 2016.04.12
  1. Correct plots:
    1. Spelling of "Instantaneous" in 1st plot
    2. Replace tau as a dummy variable
    3. Change the ordering of curves so H appears in a logical place in the list. I think it's OK.
    4. Think about reducing the complexity or changing one of the plots to be more instructive (FAIL!)
    5. Perhaps add some simpler plots with annotation such as "This curve is the difference of these two" or "This curve is the integrated difference of these two" etc. Simpler?? Lol
    6. Perhaps change the colors or the background colors to make them stand out better or have better logical cohesion (tried this: but it didn't work out too well)
    7. Think about a u related curve which represents u integrated over a window of length $T_s$ like the dashed purple $T_s$-length-sliding-window integrated $w_2$. 
    8. Annotate the x-axis better to show how a G&L period lines up with time, and at what point in the period G&L report the values (the start of the period, bizarrely).
    9. Rename those variables with subscript 1 to either drop the subscript completely (so T matches G&L's T for example), or replace 1 with "avg" for average or "mean" and likewise replace the subscript 2 with "inst" for instantaneous. Maybe on another iteration?
  2. Write out equations that demonstrate:
    1. Continuous and discrete time average and instantaneous equations, in both canonical and G&L (non-canonical) format, like in SIM8. (Also correct SIM8)
    2. Explain why you use the notation you do, and show a table with all the variable names
    3. Show how the other output variables are also affected (not just T).
    4. Write out the equations that show the solutions for all the curves.
    5. Discuss why the canonical equations make sense with a simple 1st order continuous time and discrete time pair (with a ZOH input), and show what G&L are doing differently.
    6. Either dispense with $T_s = 1$ at the start (i.e. justify dropping it), or carry it through so you don't forget that integrating over a window of length 1 is the same as a moving average over that window.
  3. Embed a spreadsheet that shows these plots and/or perhaps an animation that demonstrates exactly what's going on.
    1. Ensure that all G&L's accounting equations are satisfied where they should be, and demonstrate continuous time curves that satisfy them over any sample period, as well as discrete time tables of results.
    2. Make sure everything still works for an arbitrary ZOH input G.


  1. These charts reflect an enormous amount of very careful work! Very impressive. I believe we are both thinking along very similar lines but with a mysterious unrecognized (at least by me) factor twisted/woven in. The mysterious factor prevents me from forming a comprehensive mechanical model without jumps. Someplace I think there is a jump that I can't put my finger on.

    You write "2. The concept of finding average rates vs instantaneous rates. G, Y, T, Yd and C are all rates (or can be thought of that way). H is the only one that's different: it's a quantity." This has me scratching my head also.

    I have some thoughts about the solution but offering them now would only be intellectual noise. I'll let you know if I have something worth sharing.

    Did you do all of this in Excel?

    1. "Did you do all of this in Excel?"

      Yes. Perhaps not the best choice. I could share the file at some point too. It'd kind of nice to be able to jump from tab to tab to see the curves that change.

    2. One of the things that should stand out is that formulated as a canonical system, it all makes sense (in continuous time): the system stays at zero until the input u starts to stimulate it, at which point all the curves start to change immediately. In G&L's non canonical form, H starts to move a full one sample period prior to the input being applied! (maybe that's your "jump" you're worried about... but I don't think so, because this is only a "problem" in continuous time). Now the reason it's a "problem" at all is because of the way I've chosen to represent it in continuous time... literally the simplest and most straightforward way I could think of... but there is not a unique representation. Still, it's a very simple system, so I was hoping the most obvious choice would still make sense... however, I confess the non-canonical nature of their equations bothered me from the get-go... I ignored my concerns and low and behold it mostly works anyway (a little trouble getting the accounting equations to be satisfied when changing sample-periods in my previous work (SIM1 through SIM8)... and "sample period invariance" ... but not quite: H is rock solid, but the other curves move a bit depending on which method of change you use... or they cross a little one way or the other (output C lags Y or Y lags C slightly for example). Same goes for choosing parameters to hit steady state gain and time constant targets... sometimes a few goofy things happen... I don't think my solutions are unique: I mostly concentrated on doing exactly what I wanted to do with H and then just worrying about the steady state behavior of the others, since I knew their dynamics would be mostly driven by H (other than a jump at the start provided by G and the size of the D matrix relative to the C matrix)... so I didn't sweat the details.

      What always bothered me the most was T the total taxes collected each period (or the average rate that period)... it was simple because it basically had just one accounting equation to satisfy: G - T should be the change in H each period... so one way to construct T was from H and G directly, since G is very simple if it stays a scaled unit step function (the G&L default), and H is the easiest of the others to write down... but when you do that you get what seems like a slightly different expression for it! It satisfies that one accounting equation, but it bothered me that it wasn't the T I get from the other formulation... this exercise has let me see why that is! It's because when I formed T that way I assumed for G an associated integrated G function which was a ramp starting where H starts to become non-zero... but that's NOT what happens because of the non-canonical way G&L chose to write these equations down: instead integrated G starts ramping up (in continuous time) one sample period later than H starts to change!... See I wanted to be able to satisfy the equations at ALL times, but it doesn't work out that well with G as given. It works out great with the associated canonical u though! So now I think of G as a delayed version of U... G, T and H only satisfy the expression at the sample times whereas the u, a sliding one period wide integral of w2 and H satisfy it over any length interval offset by any amount (which is more satisfying!). I know Kirckhoff's laws (for example) in linear circuit theory work over ANY interval of time (e.g. the number of electrons that have flowed into a junction always equal the number that have flowed out) , so I KNEW there must be an expression buried in there which did so as well! (Even now I have to stop and scratch my head a little to make sure I'm writing the correct thing... I'll turn MathJax on so I can express myself more clearly).

    3. Something I wrote above make me think I still haven't plotted these curves in the best way: See the continuous time analog of the discrete time T1 curve is the w2 curve integrated over a single sample period Ts. T1 is on "G&L's schedule" as satisfies their H, G & T accounting equation at each sample time... perfect for discrete time. In continuous time that equation is satisfied over any interval of time by H, integrated-from-negative-infinity-u (playing the role of G) and integrated-from-negative-infinity w2 (playing the role of T). Can we justify those roles? Sure: as the sample period decreases to 0 the world is made right again.

      Another thing that really bothers me about G&L is that they start with what they call "Period 1" which is when all values are at their initial rest state of 0. However Period 1 refers to an interval of time, and squaring that with the usual continuous time practice of starting all activity at t=0 means that Period 1 extends from t=0 to t=Ts, but the numbers in their table give the state of the system at the BEGINNING of that interval. On my first plot (the one labeled "Canonical") I write a couple of continuous time values underneath the Period numbers on the x-axis so you can see what the underlying continuous time time-line would look like.

      So in terms of the values on the x-axis, Period 1 extends from 1 to 2, but it reports on what's happening at 1. But 1 corresponds to t=0 and 2 corresponds to t=1. Make sense? It's a minor thing I know... I could have just ignored their goofy way of doing it, but I wanted the values to match their table.

    4. What I'm trying to say (regarding accounting and continuous vs discrete time) is that if you start with the top plot ("Canonical") then the black H curve is the difference between the dark blue curve (integrated u) and the light orange (almost yellow) looking curve called integrated w2.

      In discrete time, you have to drop down to the next plot ("Non-cannonical") and then H, G and T1 satisfy it at every sample time. What's the continuous time function that joins these points for this particular sample period? That would be the dashed light-blue u and the dashed purple Ts-window-integrated w2 in the "Canonical" plot at the top.

    5. The bottom curve shows this is true (if you ignore all the noise from there being so many curves!)... the dashed purple curve falls right on top of the light-orange curve up until Period Number = 2, and then it bends and follows the solid red T1 curve.

  2. I understand what you are saying, mostly.

    Here is a quick table to express the problem in a different way:

    & G_{-1} = H_0 + T_0 \\
    & G_0 = H_0 + T_0 \\
    & G_0 = H_{+1} + T_{+1} \\
    & G_0 = H_{+2} + T_{+2} \\
    & G_0 = ............

    where G is a POSITIVE change in Government spending.

    Hope Mathjax works the first time.

    1. Sorry, Roger, I haven't turned on MathJax yet! Lol... I'll do that now. Surprisingly, you may have to reenter your table to have it show: but here's a trick. BTW here's a good MathJax trick: if you past a MathJax coded comment w/o saving it in a text file somewhere first, but you want the LaTeX code back... just press the button to delete your comment... when it takes you to the page to ask if you really want to do that, it shows you comment in it's encoded form. You can cancel the delete at that point if you want. Of course you can always "view source" but then you have to search for you stuff in a sea of other HTML.

    2. OK, great, you didn't have to redo the table! Hmmm... I'm not sure I agree though. G&L's equation holds with my $T_1$ (where the 1 subscript is not a time index.... I think I'll change that to avoid confusion), basically this: $\Delta H[n] = H[n] - H[n-1] = G[n] - T_1[n]$. Now of course G stays a constant after it steps up, so you could replace that with a constant value like you have it. In continuous time I have that written as $H(t) \equiv \eta(t) = \int_{-\infty}^{t} u(t') dt' - \int_{-\infty}^{t} w_2(t') dt'$. If you want that in terms of a single $T_s$ sized sample period difference, then $\Delta_{T_s} H(t) \equiv H(t) - H(t-T_s) = u(t) - \int_{t-T_s}^{t} w_2(t') dt'$ where $\int_{t-T_s}^{t} w_2(t')$ appears on the plot as the dashed purple curve labeled $\Delta \int w_2(\tau) d\tau$ with $\tau$ of course being a terrible choice for a dummy argument since I use various $\tau$ related names as variables I've plotted... shoot!... something else to fix.

    3. Keep in mind that some of the simplification in my equations rely on $T_s = 1$. I usually carry $T_s$ along everywhere, but not necessarily in all of the above. $T_s = 1$ means sliding integrals of length $T_s$ become the same as sliding averages (for example).

    4. About those variables with names based on $\tau$: they are what relates the average system of $T_1$ to the instantaneous system of $T_2$: and yes that really is both a bar and a dot on one: $T_1[n] \equiv \bar{\dot{\tau}}(n T_s)$ where $\bar{\dot{\tau}}(t)$ is the solid red curve. It's the continuous time average rate of tax collection over a window of time of size $T_s$. An "average rate" means it gets a bar on top of a dot. ;) (which you see Excel sometimes struggles with)

  3. $$


    & G_{-1} = H_0 + T_0 \\

    & G_0 = H_0 + T_0 \\

    & G_0 = H_{+1} + T_{+1} \\

    & G_0 = H_{+2} + T_{+2} \\

    & G_0 =



    Thanks for turning it on, but it took a few tries before my test page worked also. Still don't know exactly what the problem was.

    My Eq. 10 assumes that $ G_0 $ transforms to $H_0 + T_0$ and then $\frac{T_0}{ \theta} = Y_0$.

    Your work with tax is only different by the tax rate from being Y. We have to be very close together. I think my ignorance with notation and the reversal in control engineering time representation are the remaining gaps.

    I am trying another picture to replace the hyperbola and the idea of Object GDP. Maybe a different presentation will help. It helps organize my thinking!

    {Sometimes I must look at something 4 times to see it for the first time! LOL]

    April 10, 2016 at 2:27 PM

  4. I added a Figure 0 to the page at

    I see you did not understand my references to the Z-plane. Well, I must not have the correct understanding, so I am making inappropriate analogies. I will keep trying.

  5. How would I change Figure 0 to make it show your starting points? Or does time flow in the wrong direction?

    I drew this thinking that government does not just make a lump jump. Usually, the jump is spread over a time period. I think of "first measurement" as the first calculated Y or T.

    I think you always assume that initial H is zero, as I do. Eq. 10 will accept an initial H if it is known.

    1. Hi Roger, I only show the above starting points to be consistent with G&L. On the top figure I show what I'm assuming time is underneath the 1st and last value of "Period #" on the x-axis. I'm with you: I prefer to start things off at t=0 if possible!

      Yes, I assume $H_0 = 0$ (or $H_1$ if you're indexing by "period" instead of time, as G&L do) most of the time, but you're absolutely right, there's no problem with assuming that $H_0 \neq 0$. That's called the "initial state" or "boundary condition." I assume $H_0 = 0$ most of the time purely for convenience.

      Fun fact: instead of having $H_0 = H_{init} \neq 0$ for a simple single state system like this with a single input, you could instead add an appropriately scaled "Dirac delta function" (written as $\delta (t)$) to your input $G$. The Dirac delta function is drawn as an arrow pointed upwards on the time-line at t=0 and it represents an infinitely high spike with 0 width, which nonetheless has unit area (look it up in Wikipedia: they have a good animation). You arrive at it by taking any curve with unit area under it (like a rectangular pulse for example, 1 unit high and 1 unit wide, centered at t=0) and then taking the limit of the process of reducing its width by a factor of A while simultaneously increasing its height by a factor of A so that it maintains unit area. A system's "impulse response" (which completely characterizes an LTI system) is thought of as being in response to such an input. When such a pulse is sent through an integrator, it produces an immediate jump up to a value of the system's state (and outputs too usually)... and if the system is a pure integrator (transfer function 1/s) then it'll change such a spike into a unit step function on output... otherwise the step exponentially decays (if the system is stable).

      In discrete time the analog is the Kronecker delta function: $\delta_n = 1$ if $n=0$ and $0$ for $n \neq 0$. You can think of any discrete time system input function as being the sum of a bunch of time shifted and scaled delta functions... which with an LTI system produces as an output the superposition of a corresponding set of time shifted and scaled system impulse responses.

  6. The third reading of your material prompts this question.

    Tom, are you trying to begin from decay zero and then integrate back in time to find the larger past value?

    I am doing just the opposite. I begin at the theoretical limit and then share that limit between two owners at the first measurement time. The rates provide the time scale.

    I hope that my Figure 0 helps to understand my approach.

    1. "Tom, are you trying to begin from decay zero"... hmmm, I'm not sure what you mean. What mostly inspired this post was the last bit of SIM8 where I posted a bunch of relationships with integrals in them (integrating over a sample period $T_s$, mostly). After writing all that out, I realized it wasn't exactly true because the continuous time state space representations I presented in the upper part of SIM8 produce outputs which can be thought of as the average rate (in dollars per period) of each output, whereas the integral equations at the bottom of SIM8 only apply to outputs which are instantaneous rates at any one moment in time. It them occurred to me that I didn't know how to write the system matrices (in particular the measurement matrices C and D) for a system which produced instantaneous rates as outputs for SIM. Thus I figured out how to do that. I have yet to show the math, but I thought I'd share the results from the spreadsheet I used to check my answers. Make sense?

    2. "can be thought of as the average rate (in dollars per period) of each output" ... to be clear, that's the average over the sample period.

  7. OK, I can now answer my own question. You are not beginning at decay zero. You are beginning with the first jump.

    Now I see what looks to me like an initial start inconsistency.

    Assume that government begins with spending. For instance, to an employee. In am thinking of an income tax.

    The government gives the employee USD 20 and IMMEDIATELY takes USD 4 back in tax. The employee gets USD 16 to spend. The employee's wealth has increased by USD16.

    I think that GDP has gone up by USD20 which is why I was using GDP as my limit.

    So I think that your very first points (first measurement) for the Y, H and T curves should begin (discrete time plot) at 20, 16, and 4.

    Does that make sense?

    1. Yes, I think that makes sense. And that may very well be what G&L specify, however my goal here has perhaps deviated a bit from G&L's concept: What I'm really doing here is based on assuming a smooth continuous time system. I want to show what such a system is like and that it can in fact obey G&L's accounting equations over any interval of time (not just the given sample periods).

    2. So in particular, in continuous time, I'm solving for a case where G&L's default government spending profile (a 20 dollar unit step up from 0 dollars at t=0, and then 20 dollars per period for all subsequent periods) is a constant rate for all time t >=0, even within each period. Thus the total gross government spending integrated from $t = -\infty$ is a ramp (straight line), starting at $t=0$ at 0 dollars and then climbing at 20 dollars per period for all time (that's the dark blue dashed line labeled $\int_{-\infty}^{t} u(t') dt'$ in my 1st plot labeled "Canonical.") You can solve the equations like that, which was my goal.

  8. Hey Roger, I take it you live in Washington State... so you're on my same time. It's 9:37 PM right now, 2016.04.10, and I think I see you looking at my page on my little globe in the upper right corner. It shows somebody from that area is looking right now. I live in Goleta, CA, so I show up too.

    1. Do you live in Ellensburg, Washington? When I click on the globe, it brings that up as the origin of the viewer who's NOT from Goleta. (I didn't even know I could do that).

    2. "On continuous time analysis, it seems to me that the base assumption such as rates would apply at smallest level but we would allow the time scale to change."

      I'm not sure what you mean.

      "U(t) is the time scale? Right?"

      No, $u(t)$ is the name of the input signal to a canonical system (see my 1st plot). If a different signal is used as input (for example, offset in time), then it's nice to rewrite that different signal in terms of a u(t) with a "canonical" time basis to it if possible.

      It's easier to describe in discrete time: in discrete time SIM looks like:

      H[n+1] = A*H[n] + B*G[n+1]

      So to make that match up with a canonical system:

      x[n+1] = A*x[n] + B*u[n]

      I define a new input for SIM called u[n] such that:

      G[n] = u[n-1]

      so basically, G[n] is a delayed version of u[n] by 1-sample period: In Z-transform terms $G(z) = z^{-1} u(z)$.

      "If yes, then the initial slope should remain unchanged. What will change would be the amount shared with the second wealth holder due to the shorter time of accumulation. Right?"

      Well, u(t) really just describes the government spending function, which doesn't care what happens to the wealth. It just jumps up to 20 dollars per period, starting at t=0, and continues forever. It's gross spending, so taxes don't matter. It's integral is thus a ramp, climbing at 20 dollars per period.

      What I think (but have not yet proved to myself) is that the equations I came up with are actually independent of the exact shape of the government spending function, provided that function consists of nothing but Dirac delta functions and step functions, and it only changes on sample times.

    3. Yes, I was looking using my tablet, as I awoke from a nap. Brain strain or maybe overload! LOL

      Thanks for the dashpot series. That fits perfectly with the seismometer physics which I previously studied but not yet conquered. That will be for first thing tomorrow morning when I am fresh.

      On continuous time analysis, it seems to me that the base assumption such as rates would apply at smallest level but we would allow the time scale to change. U(t) is the time scale? Right? If yes, then the initial slope should remain unchanged. What will change would be the amount shared with the second wealth holder due to the shorter time of accumulation. Right?

      As usual, I am probably missing something(s).

      OK, delete and answer my own question.

      Below the first measuring point, the problem becomes undefined. Reasoning: The rules are the rules. They apply at first measurement point. Below that point, meaning earlier in time, we have no rules. Because we applied the rules at first measurement, we have already applied all we know. There would be no knowledge prior to what we know.

      Does that make sense?

    4. "Well, u(t) really just describes the government spending function, which doesn't care what happens to the wealth."

      Coupled with my previous just published comment (deleted and self answered), I think I finally understand the z-transform. u(t) is extending the beginning point into the past. Right.

      If so, we are extending the rules into the past. before the first measuring point. BUT, wouldn't this still be an exact projection?

      Hmmm. Time for bed. Will await morning before looking for answers.

    5. "Below the first measuring point"

      Does that mean at times prior to the first sample?

      Out first sample (according to G&L) shows all variables at 0 including the input. One thing about linear systems (even the unstable ones, ... at least in theory) is that if they are at zero and have a zero input, they will stay at zero forever. I say "in theory" because with an unstable system, the slightest bit of noise will cause it to leave 0 and start doing something (usually something bad). But even there, for a discrete time system, with quantized inputs (e.g. a digital filter), the noise could be below the quantization level, and then it would sit there at 0 waiting forever (or until a big enough disturbance destabilized it).

      So, I see a system starting at zero, and I'm constructing one possible way in which it could be changing between sample times. My solution is certainly not unique! It's just an easy one, because it's smooth. I could have tried assuming some crazy oscillatory spending pattern between sample points, but that would be too much work.

      In terms of a discrete time system, many different continuous time state trajectories map to the same set of discrete time samples, so there's definitely ambiguity there! In a sense I'm picking the simplest one.

    6. "BUT, wouldn't this still be an exact projection?"

      I'm not sure what that means.

      My plots in this post are a bit confusing I think. they are dense with information.... but they have not been pruned down to a size to maximize clarity. I wouldn't waste too many brain cells on this particular problem (unless you want to of course). Like I explain above the canonical vs non-canonical time base is kind of a side issue: the main thing I wanted to resolve for myself (which I still haven't gotten to) is outputs which represent averages over each sample period (average rates, since all outputs are rates) vs instantaneous ones. Again, that ties into my SIM8 post, the bottom part of which would be true for instantaneous outputs but is not correct (other than approximately) for average outputs (which the top of of SIM8 produces the system matrices for). So don't believe all those integral equations in SIM8: they're not correct!

      Anyway, I figured out how to translate instantaneous and average outputs, which is kind of cool. The system will require a different C and D matrix for those equations at the bottom of SIM8 to be exactly correct. But ALL these continuous time "equivalents" make assumptions due to the ambiguity in mapping a continuous time system to a discrete time system at just the sample points.

    7. "(which I still haven't gotten to)" ... gotten to posting that is. I have 8 pages of notes... Lol... but "I pity the fool" that tries to decipher them. I need to think about how to boil the information down into a clear and concise explanation. Ideally it will be half a page. It's not rocket science... but it was something that's been bothering me for a while now. I'm happy I finally got a handle on it! One of the many things about this seemingly simple 1st order system which makes it interesting. (I think some of this ties into to Jason's posts as well)

    8. Z-transforms: recall that Z-transforms are a coordinate transformation, but one (like it's cousins Laplace and Fourier) with some interesting properties.

      One of them is that a shift in one domain, is a phase change in the other. Thus if a waveform is shifted in time, it does not shift in frequency: instead the values in frequency undergo a change of phase (which can be a function of frequency). The envelope (magnitude) of the frequency function doesn't change.

      The same thing is true in reverse: a shift of a signal or waveform in frequency (in the frequency domain) does not change the timing or the envelope (magnitude) of the system in the time domain: instead it changes the phase (on a sample by sample bases in discrete time).

      You can think of how this is true in frequency: take a voice signal that modulates an RF carrier (and AM signal). That's a great example of shifting a signal in frequency. The time domain's envelope doesn't change: only the phase of the signal (which rapidly changes from sample to sample, or from time to time, corresponding to the new RF frequency it's centered at).

      Take a look at this table of Z-transform properties, and see the "Time shifting" property in particular. Thus if $X$ is the Z-transform of $x$ then if we delay x by k sample periods, that's equivalent to multiplying X by $z^{-k}$. If we're evaluating X(z) on the unit circle (i.e. we're looking at the discrete time Fourier transform of x), then |z| = 1, and this will not change the magnitude of X(z) and in particular where it's magnitude is large or small around the circle. What it does is changes the phase of X(z) at every place it's non-zero.

      Recall that all complex numbers can be represented by a phase and magnitude. That's exactly what the notation exp(a + iw) represents, where both a and w are real numbers: exp(a + iw) = exp(a)exp(iw). Thus in polar coordinates, this number (call it x) is of magnitude exp(a) >= 0, and at a phase angle of w radians.

      Thus $x^n$ = $exp(a)^n exp(i n w)$ which is a number with magnitude $|x|^n$ and phase $n w$. Multiply the magnitudes and add the phases. If a number is unit length, then all you're doing is adding the phases.

  9. ""Below the first measuring point"

    Does that mean at times prior to the first sample?

    Good Morning!

    I went to bed thinking about things I had never considered before. Exactly what do we mean about time and measurement? I had my new Figure 0 to aid me, and that helped. This is almost worth a post of it's own. We'll see where it goes.

    First an overview, then back it up.

    Begin at the first point where we know the situation without question. "Example: Gov issues USD, has someone to issue to, and charges a tax." This is a mechanical system. It is a unit, an object, something you can hold in your hand, and something that you can be put in jail for if you fail to follow the rules.

    If we separate any of those three components, we have new system with only two components or maybe four components. A new system anyway.

    OK, we have the gov issue, three component system that we are talking about. We define it at the first measuring point, when gov issues money AND worker receives it. Two actions, one by each party. Tax is extracted AT THAT POINT IN TIME. These three components can be used to establish a line on a chart, my Figure 0. More exactly, just one vertical line on the chart with a zero point, the tax share T , the receivers share H, and the total transaction Y.

    We have not yet built the whole Figure 0. To do this we need to add another rule. This becomes another system, because we now have four components. The new rule is that the receiver can transfer the USD if a tax is paid at each transfer. The receiver becomes a decision maker and can spend in any time period moving forward in time. In this four component system, if gov. issues the identical amount each period, we build this Figure 0 where we show TWO initial replicated-issues traced to stability.

    Now, what happens if we decide to begin an analysis at first issue (the single vertical line) and look BACK in time? A projection that we can do with math and imagination, but impossible in the physical world. The curve we would develop would find a limit! The limit would be at zero in the time scale. The limit would still be a vertical line. The intersection of the curve and vertical line would never occur. And the system could still be present at astonishingly small scale.

    Just to keep the time line clear, the system issue is at time 1, time zero is the start of the period preceding the issue and time two is at the end of the period following the issue.

    Yes, that concept came into my head last night as I went to sleep. Does it make sense or is it just part of a nightmare? LOL

    Oh, a surprise! The curve traced on both sides of the issue line is a hyperbola! (I think)

    1. I understand the words you're writing, even some of the phrases, and perhaps a sentence or two, but I'll wait for your post (or page) on this issue, because I'm not seeing what you're getting at. And I don't see where a hyperbola comes into play. With the default input, every curve produced by SIM is (a scaled version) of one of these types (only considering t >=0, with all assumed to be 0 for t<0):

      1) 1 : a unit step

      2) t : a unit ramp

      3) 1-exp(-at) : a unit exponentially delayed rise, where parameter "a" is determined by the system's pole

      4) A linear combination of curves 1) and 3)

      5) Extending to my set of curves, then a linear combination of curves 2) and 3) as well (for example the integral of w2)

      Now you could take a product or ratio of any two of those, but I'm still not seeing a hyperbola in the mix.

    2. ... It's actually simpler than what I wrote above: With the default system input, SIM only produces two kinds of outputs:

      1) A unit step (call it U(t) to distinguish it from general input u(t)). This is because this is the form of the input function cast in canonical form and moved to continuous time.

      2) exp(At)*U(t): An exponential times U(t). The complex frequency of the exponential is A which is determined completely by alpha1, alpha2 and theta (the parameters). In this case it's a 1st order system so A is a scalar and the negative of it's reciprocal is the system's time constant (Tc) (i.e. it's decay time). There are NO OTHER "frequencies" in the system. A happens to be negative, so this exponential is decaying.

      All SIM curves (in the canonical, continuous time version of it that I've picked as in some sense the simplest version) are linear combinations of the above, meaning p*U(t) + r*U(t)*exp(A*t) where p and r are scalars.

      Moving on to a discrete time version of this, the outputs are ONLY samples of these continuous time curves. Moving on to SIM's non-canonical setup: they are potentially delayed by one time step as well (H isn't, but the others are).

      All other curves based on these that I've presented (time derivatives or time integrals) at best only introduce higher orders of the step function (t, t^2, t^3, etc) and more decaying exponentials: U(t)*exp(At). I only integrate once at most, so that gives me a ramp (t), but those curves are derivative curves I used to help me understand how to make instantaneous relate to average.

      Making these statements more general, if your input is in general u(t) (not necessarily just a step) then this system will produce as output ONLY linear combinations of:

      i) u(t) itself

      ii) U(t)*exp(At) (yes, that's big U, the unit step function)

      The only difference is that U(t)*exp(At) (AKA the system's "impulse response") could be delayed as well as scaled.

      Then the system's outputs including its state (H) will be linear combinations of scaled and delayed versions of ii) and scaled versions of i).

    3. ... and Good morning Roger!... (sorry I jumped right in w/o a salutation ;)

    4. The hyperbola is almost a distraction. It might be useful in the future, if we can link it to the vertical start line, a zero time and the negative (which would be when the government replicated period-issue stops). But a hyperbolic distraction for now. LOL

      I think I see what you are doing. You are continuously scaling time to smaller and smaller increments. (or larger and larger increments). Then each rate of change rule continues to apply (i.e., the slope between time periods is the same) so the magnitude goes down (or up) and the location on the paper projection shows the proper curve. The same as my Figure 0. Right?

      The hyperbola relates to my Figure 0. Each of the curves is part of a hyperbola (I think). The formula is $ k = \frac{T}{Y} = T ( \frac{1}{Y}) $ where k is a constant and the last term is the more familiar appearance of the hyperbolic formula 1 = a*b.

    5. and term "T" could be time or taxes because they are proportionate under our rule of "tax per period".

    6. I am forming an opinion. You and I are doing the same thing in the sense that we are building the same SIM curve. You are working on changing the time line but holding the slopes the same to find the magnitude at the second (and later) measuring points. I am changing the slope between measuring points.

      Since I am changing the slope of the linear line between the two measuring points, I need to hold the time increments constant but change the rate-of-change. To do that, I need the H term. My formula gets a new H for each period iteration so we effectively iterate column by column (in the spread sheet) using term "G + H" where "G + H" becomes a constant over several columns of iteration. (but we do not need to turn the spreadsheet iteration feature).

      Does this all make sense now?

    7. Hi Roger, the expression T/Y as a function of time will be of the form:
      $$\frac{T(t)}{Y(t)} = \frac{p_T + r_T exp(a t)}{p_Y + r_Y exp(a t)}$$
      where $p_T, r_T, p_Y, r_Y$ and $a$ are all constants. I don't think that works out to be a hyperbola. But if you can show it, I'd love to see. Perhaps another post? Or perhaps you're speaking of their steady state values $\bar{T}$ and $\bar{Y}$ as a function of one of one or more of the parameters? Equations (4) and (5) of my SIM7 show that in steady state we have:
      $$\frac{\bar{T}}{\bar{Y}} = \theta$$

    8. ... and since I know ahead of time that $a < 0$ (i.e. the system is stable), then I can tell that $p_T / p_Y = \theta$ as well.

    9. ... so your $k = \theta$ in steady state.

    10. ... but I think you already knew that.

    11. You write " With the default input, every curve produced by SIM is (a scaled version) of one of these types (only considering t >=0, with all assumed to be 0 for t<0):

      1) 1 : a unit step

      2) t : a unit ramp

      3) 1-exp(-at) : a unit exponentially delayed rise, where parameter "a" is determined by the system's pole

      4) A linear combination of curves 1) and 3)

      5) Extending to my set of curves, then a linear combination of curves 2) and 3) as well (for example the integral of w2)"

      I would write for my approach:

      1) 1 : a unit time step
      2) H ; a memory of distance already traveled up the vertical magnitude scale
      3) use the rates of change to find the new magnitude of Y using (G + H)
      4) Y : record Y in the spread sheet. Use Y to calculate a new H
      5) repeat steps 1 through 5 using spreadsheet columns until stability is reached.

      Hmmm, Am I still correct in saying that I am working on the magnitude axis and you are working on the time axis?

    12. Roger, I'd use my later comment where I boil that down to just two things:

      1. $f_1(t) = U(t)$ (a unit step)
      2. $f_2(t) = U(t) e^{a t}$

      and linear combinations of 1. and 2. (i.e. $p f_1(t) + r f_2(t)$ for any scalars $p$ and $r$).

      And yes, I'm thinking of functions in time, but I'll switch to a particular time if you like, or to steady state, or even to the frequency domain. ;)

    13. Recall that $e^{a t}$ is a very interesting function: you can differentiate it, integrate it, or shift it in time ($e^{a (t + t_0)} = e^{a t_0} e^{a t} = k e^{a t}$) and all you get as a result are scaled versions of $e^{a t}$ back and maybe an added constant (if you're integrating starting from some fixed point in time up to time $t$).

    14. ... assuming you differentiate and/or integrate wrt $t$ of course (and that $a$ is a constant).

    15. In my 9:40 AM comment above, I forgot to multiply the first expression of T(t)/Y(t) by the unit step U(t).

    16. "2) H ; a memory of distance already traveled up the vertical magnitude scale" I'd put that by saying H is measured in dollars without respect to a sample period. That's why the y-axis of my three plots in the post say "dollars or dollars per sample period"... the plain "dollars" part of that is for H, but since I integrate several of the variables and plot them as well, it also applies to them (e.g. integrated u and integrated w2).

      "3) use the rates of change to find the new magnitude of Y using (G + H)"

      I agree.

      "4) Y : record Y in the spread sheet. Use Y to calculate a new H"

      I disagree with that. H is not dependent on Y at all. Y factors out of the equation. H only depends on $G, \alpha_1, \alpha_2$, and $\theta$.

      "repeat steps 1 through 5 using spreadsheet columns until stability is reached."

      Other than the final step number you give, I agree with that.

    17. ... I guess you could say that Y is hidden in those parameters and G and previous values of H. But the same could be said for T, YD or C. I think of all four of them as outputs (linear combinations of H[n-1] and G[n]), not contributing to the dynamics of the system at all... simply measurements of it. That's maybe not how reality is, but that's how this model is. G is the engine, and the model of H (dependent on the parameters) is the system reacting to it. The others (Y, T, YD and C) simply form linear combinations of those two and report on them. There's no feedback path from them into the rest of the system. That's why all curves have $e^{a t}$ in them: because the state update model of H does. You won't find any $e^{b t}$ in there where $b \neq a$ because there ARE NO other characteristic frequencies, because the model for H has no other characteristic frequencies. The only way you could get a $e^{b t}$ in the outputs is if the exogenous input G had that form.

    18. "Hi Roger, the expression T/Y as a function of time will be of the form:"

      Hmmm, My T needs the included term $ \theta $ . See if I can explain:

      Use the example of Gov pays employee (Y), collects tax (T) using tax rate $ \theta $ per year. This is where the system becomes important.

      Y is an annual measure of output so can be considered a rate in the macro economy.
      T is an annual actual extraction of wealth so can be considered a rate in the macro economy
      $ \theta $ is a ratio of tax to output so can be considered as a rate per rate in the macro economy.

      So, when you write $ \frac{T}{Y} $ you really mean $ \theta $ .

      And when I write Y/T is a hyperbola, I really mean $ \frac{ T }{ \theta } $ is a hyperbola.

      The tax rate becomes a rate per rate on the time scale.

      Does this any of this make sense?

    19. ... at it's heart, a state space model is a model of the states. The state update equation (not the measurement equation (i.e. the y equation)) is the most important part. That's why the matrix A is of utmost importance. The measurements provide no feedback to the system, they simply form linear combinations of the states (x) and the exongenous input (u) and output them. The only way those outputs affect the states of the system is if you wrap a feedback path around it to create a new system (which is what feed back control engineers do). In SIM the only state is H. There's no way to even transform the system into one in which the states are anything other than H or a multiple of it, and that's because you can only transform a system's state representation with linear combinations of it's current states. So if the current states are represented by the Nx1 vector x, then you can make an equivalent system with states V*x where V is some NxN non-singular constant matrix. But since x is a scalar in SIM that means V must be too. So only H or a straight multiple of it can represent the system state of SIM.

    20. "θ is a ratio of tax to output so can be considered as a rate per rate in the macro economy."

      Yes, but unlike T or Y it's a constant system parameter. T/Y = θ, yes, I agree that's true, even though both T and Y have the form (as a function in time) of U(t)*(a constant plus an exponential), they both are a scaled version of this same expression (I'd forgotten that this was true, but looking at my SIM7 confirms it). So yes, I agree that T/Y = θ for all t.

      "And when I write Y/T is a hyperbola, I really mean T/θ is a hyperbola."

      Hmm, I think you had it correct before. If you take w = T and v = 1/Y and plot w vs v it will be a "rectangular hyperbola" (as Nick Rowe calls it) since w*v = θ = a constant.

    21. ... it's "rectangular" because its asymptotes are the x and y axis themselves (rather than lines that cross at some other angle ... other than 90 degrees that is).

    22. OK Roger, I have to go now, but I'll catch up with what you've come up with later. Have fun!

  10. You write "G is the engine, and the model of H (dependent on the parameters) is the system reacting to it. "

    Yes, that is conclusion I am coming to.

    1. This comment has been removed by the author.

    2. Great! Like I wrote above, the matrix A is the heart of the system. If A is dimension NxN (meaning the state vector is dimension N), then there are at most N distinct characteristic frequencies of the system: $\{a_1,\, a_2,\, ... \,,a_N\}$ (where the $a_k$ are (possibly complex) scalar constants) which produces corresponding outputs of the form $u(t) e^{a_1 t},\, u(t) e^{a_2},\, ... \,,u(t) e^{a_N t}$ all summed together. I.e. the output will have the form:
      $$y(t) = u(t) \sum_{n=1}^{M} c_n e^{a_n t}$$
      for some set of constants $c_n$ and $M \le N$. Now that I've said that I have to tell you it's not actually true, because the system could have other outputs that look like $t^m e^{a t}$ added in, where $m \le N$ but only in "degenerate" cases. Such a matrix $A$ producing those cases is a tiny perturbation $\epsilon$ away from another matrix $A'$ which does not produce those degenerate solutions.

    3. ... I should say the output y will have the form of that sum for an impulse response (or for an initial state vector (same thing)), and assuming the $D$ matrix is $0$. If $D \neq 0$ then you have to add in a scaled version of input signal $u$ as well to $y$.

      Think of a linear system as a collection of tuning forks, all welded together in parallel and tuned to different frequencies in such a way that when you hit it against the edge of a table you get its impulse response (all those frequencies mixed together), and if you hit it a bunch of times in rapid succession with different force each time you get a bunch of frequencies added together and a bunch of impulse responses from different times added together. You can approximate any input $u$ as a bunch of different responses to different strength bangs of this compound tuning fork against a table at different times. (A tuning fork doesn't stop ringing in response to a first bang against a table just because you banged it against the table a second time and third time).

      SIM is a single tuning fork, and the atmosphere is so viscous that it doesn't really ring when you bang it, it just produces a decaying exponential.

    4. "You can approximate any input u as a bunch of different responses to different strength bangs..."

      I should have written:

      "You can approximate any input u as a bunch of different bangs (against the table) of different force at different times."

      And really I'm describing a linear time-invariant (LTI) system. More generally the analogy holds, but if the linear system is time varying, it means the tuning forks are actively being tuned to different frequencies as functions of time as you bang them against the table. One complicated compound tuning fork!

    5. ... that approximation of u (as a series of different strength bangs) becomes exact as the time between bangs goes to zero, and the force of each bang becomes vanishingly small.

    6. ... an unstable system is one in which the oscillations (at least one of them) of the compound tuning fork grows in time rather than diminishes... so eventually it breaks apart!

    7. ...that $u(t)$ above, multiplying the summation, is supposed to be the unit step. Sorry I keep changing that! It's there just so it's clear that the tuning forks don't start ringing before you bang them the 1st time.

    8. ... if you want to know what the characteristic frequencies are of any random square matrix A, then type this in Octave at the prompt:

      A = randn(5)

      which produces a random 5x5 matrix and displays it on the screen. Then type:


      And it will list A's characteristic frequencies = poles = eigenvalues. If you want to see what happens with a degenerate A, then try this:

      A = [0 1; 0 0]

      And you'll see it print on the screen:

      A =
      0 1
      0 0

      Then try the eig() function again:


      eig() should fail on matrix A because it's degenerate (i.e. it's not "diagonalizable"). In general eig() will fail on any matrix A which is composed of Jordan blocks (at least one) of dimension > 1, or any derivative of such a matrix of the form B*A*inv(B) where B is the same dimension as A and invertible. However if you construct such a matrix, then you can add a tiny perturbation to it, and it will be diagonalizable again:
      A = A + 1e-16*rand(N)
      will be diagonalizable with probability = 1

      Sadly, degenerate cases are not as rare as you'd think because they derive naturally from a system with a state vector composed of position and it's derivatives (velocity, acceleration, etc) in a frictionless (or air resistance free) environment. However if you were to try to find one by checking random matrices like this: eig(rand(N)) (for any positive integer N), you would probably never find one. So they're rare in the sense that out of the pool of all possible matrices, they're a tiny minority.

    9. "And it will list A's characteristic frequencies = poles = eigenvalues."

      As you may know by know, I could add to that list of synonyms:

      "And it will list A's characteristic frequencies = poles = eigenvalues = roots of its characteristic polynomial = solutions s to {det(A -s*I) = 0} = list of amplifying factors for every vector (v) which is only scaled by A (i.e. vectors (eigenvectors of A) v such that A*v = a*v, for a scalar a (where a is the eigenvalue))."

    10. ... a couple of more to add to the list of synonyms:

      1. modes
      2. resonant frequencies

      It's kind of strange thinking of a system like SIM having a "characteristic frequency" because it doesn't oscillate (like a tuning fork does), but that's only because the frequency is all real (no imaginary component). If a frequency has an imaginary component it will oscillate. The bigger the magnitude of the imaginary component, the higher the frequency of oscillation it will have. The real part controls the shape of the envelope of those oscillations... for stable systems that's a decaying exponential. So as the magnitude of this real part gets bigger, so too does the rate of decay of this exponential. So in general, the further the frequency is from the origin (i.e. the overall magnitude of the complex frequency) the "faster" the system will be: it'll either decay (or rise) fast or it'll oscillate fast or some combination of the two. Conversely, the closer a frequency is to the origin (on the complex S-plane) the slower it'll be. Rules in the Z-plane are similar but a little different, but I'll talk about those some other time.

    11. ... also, the only reason (from a system theory perspective) a system like a tuning fork seems to produce purely real oscillations, is that it actually has two conjugate "natural frequencies" which respond in such a way that the imaginary components cancel each other out. Thus for every real valued oscillation a system produces it needs two extra dimensions of its state vector. The smallest dimension system producing strictly real valued oscillations you can have is thus a 2-dimensional system having a 2x2 A matrix.

    12. Well, I just tried this myself:

      A =
      0 1
      0 0

      Then try the eig() function again:


      and it didn't cause an error, but if you try to find the diagolizing vectors V like this:

      [V,D] = eig(A)

      and then try to diagonalize A with those:


      you'll notice that it doesn't work (and I get a warning too): the result is not a diagonal matrix

    13. then go back and try that with any random matrix, and you'll see that V does diagonalize A. For example, for a random 3x3:

      A = randn(3); [V,D] = eig(A); inv(V)*A*V

      Just cut and paste that and you should see a diagonal 3x3 as a result, which if you now type


      You'll see is equal to D. I.e. D == inv(V)*A*V. The eigenvalues of A are on the diagonal. You can't do that with degenerate matrices. Their V matrices are not invertible, because not all the columns are linearly independent.

  11. Tom, I liked your "dashpot" explanation. It really helped a lot. I am still integrating it because there is a lot of content here.

    I think I see why you have a problem with the idea that we can integrate the Y term. The example dashpot has a delayed force in relationship to the distance traveled. An integration would be a force ahead in time or at least instant, as if the actuator did not directly push on the dashpot.

    Here is another dashpot analogy that I think applies:

    If this word picture doesn't work, I can easily post a picture. Hmm, I had just better start with a picture. That will take a little longer. The idea will be that the dashpot is not directly connected to the actuator. Instead, the actuator is connected to a pivoted lever. One end of the lever is hinged, the other end is connected to the dashpot mounted against the unmovable wall as before. Now, any motion of the actuator will be amplified when it travels to the dashpot. Place a distance scale under the dashpot and a second scale under the actuator. One unit of actuator motion will cause K*dashpot motion. The force will be inversely related.

    I think I need to say something like that in my state-space representation.

    Thanks so very much for developing the dashpot analogy. It is great!

    1. You're welcome. I eagerly await your diagram.

    2. Also I added a tuning fork analogy for you in the above comments... a compound tuning fork with many (N) tuning forks mounted in parallel as an analogy for bigger dimensional systems with state vector x an N-vector, N > 1.

    3. Hi Tom, I added the dashpot diagram Figure 0.1 in the usual place at

      Next, I will think about updating my state-space diagram and how to write sensible transform equations. I am still weak here but improving.

      BTW, Brian Romanchuk has a new post on SFC models.

    4. I give up on the ridiculous hyperbolic argument. A serious overlooking of the unit differences on the vertical axis and horizontal axis. Sorry for the distraction.

    5. Thanks Roger. I left some questions and a request for you there.

    6. "I give up on the ridiculous hyperbolic argument."

      You want ridiculous??? Try my brand new Figure 2 above. The curves are the same as Figure 1, but I tried to make it even more confusing. What do you think? Did I succeed?? Lol.

    7. Well, having stared at that mess (that I so painstakingly constructed) for a while... I do have a glimmer of hope of how to **SIMPLIFY** it so it actually makes sense. Wish me luck! ;) (but that will be another day I'm afraid).

      I think I can tie the G&L variables in as integrated values over the preceding period. The *real* values are the rates. So, for example, u(t) is the real government rate of expenditure (but I'll call it something better), and G is that integrated over the previous sample period. Same for taxes, etc. I think I can keep the purple dashed curve, but I won't move it: it stays fixed like H. I'll ditch the red curve. I'll keep the green curve in just one plot. I think I can go down to just two plots... maybe even one. My new labels will be:

      "Rate of government expenditure"
      "Total amount government spent"
      "Amount government spent over preceding Ts sized interval (G)" (this will actually be a new curve)
      "Rate of tax collection"
      "Total amount of taxes collected"
      "Total amount of taxes collected over preceding Ts sized interval (T)"

      and of course:

      "Total cash in economy (H)"

      So seven curves total with **SIMPLE** labels. Yes, I think one chart should do the trick. All the government spending ones can be shades of green, and all the tax ones can be shades of red, and H should still be black.

      That way it will be logical and coherent! ... believe it or not, I set out to simplify early today, and I did exactly the opposite. It wasn't till I had a chance to look at my "handy work" with fresh eyes that it seemed funny to me.

    8. "Roger SparksApril 11, 2016 at 11:16 AM
      You write "G is the engine, and the model of H (dependent on the parameters) is the system reacting to it. "

      Yes, that is conclusion I am coming to. "

      Good Morning Tom,

      Several things:

      First, Term H as the " SYSTEM ". Maybe I wrote in error, but I think G + H is the SYSTEM . Term H by itself is MEMORY for the system but term G is also MEMORY . In my view, we could think of G and H as a husband-wife team where both have memory and decision making abilities.

      Second, I will answer (or try to answer) your question about the dashpot Figure 0.1. You are seeing it a little differently from the way I see it, so I will try to account for that. In particular, I think you focus on the "state" of the system where I focus more on the system as a dynamic entity. My analogy of the husband-wife team illustrates that focus.

      I will be less noisy today. Several other things to do, and I really want to better understand the state-space notation and matrices better to reduce the danger of intellectual noise. I am afraid I have already generated too much intellectual noise.

      Next, I will work on the dashpot on page Mathjax.

    9. Tom, I updated the dashpot diagram and added text. Would you disagree with any of this? Or maybe there are features misrepresented?

    10. Two points:

      Are there any features not shown like velocity? I think not, all the money is always observable on the money scale.

      2nd, If government also has a memory, then the system can be replicated many times with G = 0 and H = 75. A stable system. Then government can suddenly 'remember' that it has USD 20. Then the system becomes dynamic beginning with G = 20 and H = 75. A few replications later, we will have stable G = 0 and H = 95.

      Do these points make sense?

    11. This comment has been removed by the author.

    12. Hi Roger, I'll have to touch base with you later... lots to do today!

    13. Quick note on including the input as part of the system: that's fine, but then you lose an input, other than it's initial state. It's like this (in Laplace transform domain):

      Y(s) = P(s)U(s)

      You can form the product of P(s) and U(s) and call that your system if you like. Recall that the Laplace transform of a unit step is 1/s, so in effect you'd just be multiplying your system by 1/s. So we already know that X(s) = bU(s)/(s-a), so then looking at one of the Y outputs, we can write as a SISO system:

      Y(s) = c*X(s) + d*U(s) = c*b*U(s)/(s-a) + d*U(s) = (c*b/(s-a) + d)/s

      I'll let you do the math to make that a ratio of two polynomials in s.

      Now this new system has two poles: one at s=a and one at s=0. It'll be the same for all the other elements of Y. The only thing that changes is the location of the zero and the gain.

      You can do the same exercise in the Z plane, keeping in mind that a unit step there is 1/(1-Z^-1), and the system P will also be different in terms of Z.

    14. ... with linear systems, you can even put the "input" on the output side of the system: order doesn't matter. This is what Jason did. The Laplace transform of the delta function (impulse) is just 1. So you could order the system like this:

      Y(s) = U(s)*P(s)*1 = P(s)*U(s)

      And everything is the same.

    15. and yes I think your comments about memory are correct.

      Also, concerning yourself with the state is concerned with dynamics: recall that A determines the state transition matrix.

    16. ... the location of poles and zeros of a system doesn't tell you about the overall gain. Imagine P(s) = N(s)/D(s) (for numerator and denominator). Then k*P(s) has the same set of poles and zeros for any scalar k. What you can tell is the relative magnitudes associated with each pole, by factoring P(s) as follows:

      P(s) = k1/(s-a1) + k2/(s-a2) + ... + kn/(s-an)

      Then you'll know the relative magnitude associated with characteristic frequency a1 vs a2 (for example), which will just be k1/k2.

      So for each of the four outputs, you know they'll all have a pole at s=a and a pole at s=0, but you can work out their relative magnitudes (again taking the full system as P(s)*U(s), where U(s) is my government spending function in the frequency domain).

    17. BTW, specifically what Jason did in his simulation was the following: to try and avoid confusion, say the unit step is unitstep(s), and the input is U(s), then he took the 20 scaling factor off of U(s) and put that in the delta function, so he had:

      Y(s) = unitstep(s)*P(s)*20

      evaluating from right to left. But of course that produces the same output. However if you follow the logic of his code, without seeing the unitstep at the end, you'd just see:

      P(s)*20, and you'd notice that all your intermediate values at different times are different than his. It's not until he applies the final integration (accumulation) of unitstep(s) at the end, that it becomes exactly the same.

    18. Re: velocity: Well the continuous time version has $\dot{H}(t)$ in the problem statement. You can always know what $\dot{H}$ is by forming A*H + B*G. Same as for the dashpot system (and it's displacement x). In both cases the derivative of the one and only state isn't a state itself in the system. However, a second order system would likely have those as states.

    19. ... and it's not until you had a 2nd order system that the discrete time equivalent would also show an explicit velocity state appearing at all.

    20. ... why aren't they states in 1st order systems? Because it's not an independent degree of freedom in those cases.