A Novice Looks at Linear Space-state

In  A Master Equation for SIM Models Using the GDP Object DRAFT we have an equation that looks (to a novice) like a regenerative amplifier combined with negative feedback.

[4/10/2016 4:39 PM edit. Added Figure 0 and  the text above it.]

Figure 0 is an attempt to picture the time shape of measurements Y and T from successive sustained jumps in government spending. A representation of the method of finding the values at first measurement is included.


Figure 0. Two successive sustained jumps in government spending, with a sketch of first measurement calculation.

[[4/12/2016 6:38 AM edit of Figure 0.1, added text following]
4/11/2016 8:59 PM edit. Added  springs to Figure 0.1]
[4/11/2016 7:15PM edit. Added Figure 0.1 and  the text above it.]
Figure 0.1 is a crude model of an economy with two memories, two actuators, and existing between two bounds. The sectors are connected but not on a one:one basis.
Figure 0.1 A physical analogy of five observable features of a two sector economy.

Five components or relationships are illustrated here:

1. Economies have memories represented by money. As a crude description, we could say that money is a memory that someone worked for government and still has the paper to prove it. The springs store store energy in the physical world. In Figure 0.1, we are thinking of energy-in-storage as being money-in-storage.

2. Economies have decision makers. The dashpots are actuators, either of which can move the state of the system. The SYSTEM money level does not need to change to move the state of the system; only the distribution of the money level between the two storage devices is required to change the state of the system.

3. When the WALLS move, the energy level of the entire system will change. For example, when government borrows from itself, it is equivalent to moving the walls of the system. The energy level in both springs will change and the position of both actuators will change.

4. The government and private sectors have different characteristics. Changes in energy level or money supply levels have uniquely different effect on each sector.

5. The system is always at a balance point where the state can be measured and located. Between measurements, the exact positions are unknown.

[hour later update] Because both sectors have memory, the system can be stable at G = 0 and H = X for many observations. If government suddenly 'remembers' that it 'has' money, we find the wall has moved to see that G = NM, H = X. The system will be dynamic until we see stability at G = 0, H = X+NM.
_____________________________________________________
|Begin Space-State discussion_____________________________|
A novice's version of a space-state representation is in Figure 1.

[4/9/2016 3:10 PM edit. Replaced Figure 1 with a new conception]
[4/9/2016 6:50 AM edit. Replaced Figure 1 with a new conception]
[4/9/2016 3:47 AM edit. Replaced Figure 1 with a new conception]
[4/8/2016 5:02 PM edit. Replaced Figure 1 with a new conception]

Figure 1. A state-space representation of Eq. 10 in the reference.

Here is Equation 10 from the reference post

$$AGDP = \frac{ G + H_0}{α1(1-TR) + TR}$$

slightly edited. $H_0$ was $H_{-1}$. The accumulated memory of $H_{-1}$ is done in a separate spreadsheet cell, so, it would be appropriate to relabel $H_{-1}$ for purposes of the state-space representation.

The reference post is still in draft mode. The correct presentation of term H is extremely important.

[4/8/2016 5:02 PM edit. Deleted the 10:22 update. Added text. [10:22 AM update]

Hmmm. Trying to write this certainly causes one to think. I believe the 1/S box covers the entire time period between beginning the calculation and final output. Therefore, output would be time $GDP_0$, with GDP created from start based on $G_{-1}$ (but incorrectly in the spread sheet taken from $G_0$ column) and $H_{-1}$.


Later, in other rows, we update the H accumulator which is another cell. This is a time memory.

Equation 10 is certainly controllable. Use Yo as a zero value input to indicate that Y(t) is part of the input. In reality, all of the continuation of Y(t) is in H(t). G(t) is an increment to Y(t) unless T(t) equals G(t).

This is certainly an interesting puzzle!


77 comments:

  1. Hi Roger, very good! I'll have a longer look at this later... I see some potential problems, but let me think about that a bit before I discuss it with you.

    Why it says mathjax in the title: here's my guess: perhaps when you first created this page, for the page title you typed "mathjax." It doesn't much matter how long you had that as a title, blogger will use that first title you give the page to name the URL, even if you later change the page title, you're stuck with the 1st URL name for good (AFAIK). Unless, of course, you delete the page, and start over.

    A couple of problems that catch my eye right away:

    1. On the left side at the first summer junction: it appears to have two outputs. Why?

    2. On the extreme right, close to the Y output, there's a junction with no indication of what's going on there. this would be fine if an input signal were being split into two or more copies, but that junction has two inputs.

    3. The names in the boxes got messed up a bit with the location of subscripts, etc.

    I think I see what you're trying to do: you want to treat some of the items I think of as fixed parameters as inputs.

    If you instead think of them as fixed parameters (I'm specifically referring to alpha1, alpha2 and theta), then the block diagrams (discrete and continuous time) you get for SIM are these, where I could have replaced the integration block with a 1/s symbol as you have it.

    Those diagrams aren't perfect. In particular the g should really be a g "dot" (i.e. g' = dg/dt) since in that model it's a rate (as I explain in my SIM4 post) expressed as dollars per period. I actually had it like that at first, then erased the dots... oh well, I haven't gotten back to correcting it. Actually all the names should appear as I have then in SIM4 for the continuous time case.

    I'll think about what you're trying to do and get back to you. But I think what you're going to find is that you're essentially treating alpha1, alpha2 and theta as time-varying inputs (along with G (or g')). In that case you have a time-varying linear system: in other words instead of A, B, C and D you'll have A(t), B(t), C(t) and D(t).

    ReplyDelete
    Replies
    1. You'll not that the flow in my diagrams is from right to left... that's a habit control engineer get into because it better matches the "flow" of the equations as they're usually written:

      x' = Ax + Bu
      y = Cx + Du

      ... but certainly left to right is more "normal" for the average reader.

      Delete
    2. Also note that in my diagrams, the parts that produce the outputs Y, T, YD, and C are all very similar, so I've only shown the one output T. You can think of T as representing the whole bunch of them (say as a vector of outputs).

      Delete
    3. You'll also note that my diagrams also implicitly represent a time varying linear system, since I spell out in the box exactly what my gain boxes are in terms of the (potentially time varying) parameters alpha1, alpha2, and theta. I could explicitly draw an signal input line for each of those parameters leading to each box which depends on them, however the system would then not be linear (even time-varying linear) wrt those "input" signals. Do you see why? Only wrt G is the system a linear system.

      Delete
    4. This comment has been removed by the author.

      Delete
    5. Also my expression for $D_T$ is correct (in the text box in each of those diagrams), but it can be simplified as I have it here in the expression for $\dot{t}$ in SIM7. That lower case $t$ is not time, it's my continuous time version of the tax variable T. In continuous time $\dot{t}$ is the tax rate in dollars per period and $t$ is the accumulated total taxes collected. Anyway, as you can see, $D_T$ should just be:

      $$D_T = \frac{\theta}{X} = \frac{\theta}{(1-\alpha_1(1-\theta))}$$

      in both discrete and continuous time. Make sense? In SIM7 I implicitly give all the components of the measurement matrices (both are 4x1) $C$ and $D$ in simple forms (i.e. I've done the algebra to simply them there).

      Delete
    6. Wow Again. I am reading this with much more understanding! But of course, still a complete novice.

      I printed out your state-space representation but just starting to look at it.

      I changed SIM_MECH in column 11 to get average G for the period. This because the column calculation will be for the entire period, and a jump using either the start or end values would be incorrect. I think this goes to your comment about using G as a time varing signal.

      I need to think about using H and the other three as time varying inputs. They can all change on a moments notice in the real world.

      I see that the spreadsheet can be downloaded in Excel format. For the jump control, I had a nice checkbox that disappeared in Google sheets. The 0 and 1 are a work-a-round.

      Delete
    7. H isn't an input. $H_0$ is an input: $H(t=0) = H_0$.

      Delete
    8. Also note that on the discrete time block diagram, I could have labeled the block called "Delay by one $(T_s)$" with "$z^{-1}$" instead.

      Delete
    9. This comment has been removed by the author.

      Delete
    10. The system is not linear wrt $\alpha_1, \alpha2$ and $\theta$ treated as inputs, however, you could create a linear system with 5 inputs based on those three parameters and the current input $G$ as follows:

      $$u_1 = B G$$
      $$u_2 = D_Y G$$
      $$u_3 = D_T G$$
      $$u_4 = D_{Y_D} G$$
      $$u_5 = D_C G$$

      Where $B, D_Y, D_{Y_D}, D_{T}$ and $D_C$ are all functions ONLY of $\alpha_1, \alpha2$ and $\theta$, as I give in the table in the lower left of my spreadsheet here (SIM6).

      Delete
    11. ... the only problem being, what do you call inputs $u_1, ... u_5$? I'm not sure they correspond to anything very intuitive or name worthy, but perhaps they do. If you figure it out, let me know!

      Delete
    12. ... one nice thing about forming those five new u inputs is the D matrix goes away, and the B matrix becomes very simple: B = [1 0 0 0 0].

      Delete
    13. This comment has been removed by the author.

      Delete
  2. This comment has been removed by the author.

    ReplyDelete
  3. "one nice thing about forming those five new u inputs is the D matrix goes away, and the B matrix becomes very simple: B = [1 0 0 0 0]."

    Yes, that is towards what I was thinking. In term B, the 1 would be for G but should there be a 1 or $\alpha 2$ for the wealth term?

    Am I running time in reverse from you?. I an thinking we get inputs at period start. Then we need to run those parameters for the period and at period end, the output GDP is the total output by both government and private sectors. Wealth is reduced by a holding factor so only part of wealth is spent in the period. All of G is spent.

    ReplyDelete
    Replies
    1. This comment has been removed by the author.

      Delete
    2. This comment has been removed by the author.

      Delete
    3. Actually, that sentence of mine wasn't quite right: I had B correct, but D would be like this in terms of an Octave/Matlab command:

      D = [zeros(4,1) eye(4)]

      If you cut and past that line into Octave it will create that matrix and assign it to D and (if you don't add a semicolon) will display the result on the console. You should get the following:
      $$
      D = \begin{bmatrix}
      0 & 1 & 0 & 0 & 0 \\
      0 & 0 & 1 & 0 & 0 \\
      0 & 0 & 0 & 1 & 0 \\
      0 & 0 & 0 & 0 & 1 \\
      \end{bmatrix}
      $$
      So a 4x5 matrix which is a 4x4 identity matrix (I) with and extra column of zeros appended on the left side (just like the Octave command above states).

      Delete
    4. Roger, you write:

      "In term B, the 1 would be for G but should there be a 1 or α2 for the wealth term?"

      Actually, the 1 isn't for G: it's for the product $B G$, which of course is a function of $\alpha_1, \theta$ and $G$ in discrete time. For the continuous time case, it's also dependent on $alpha_2$.

      ... and now that I think about this a little more, I don't think this idea of forming these five inputs is going to work as a linear system, if you want it to look like a linear system from the perspective of something other than just G. The reason is the matrix A (a scalar in this case) is still a problem: we can't have A itself as an input because it's multiplying the state (H). That's not a linear system. The state is never multiplied by an input signal.

      Thus I don't think you'll ever have a situation where you can describe such a system with the linear state update equation which involves a state transition matrix $\Phi(t+T,t)$ which transitions the states from time $t$ to time $t+T$. That transition matrix obeys the laws of superposition, and without superposition you don't have a linear system. Now we might be able to approximate a linear system about some nominal state solution and input signal (by using the usual technique of calculating a matrices of partial derivatives (i.e. the "Jacobians"), but I just don't think there's a way to write this as a truly linear system except with a single input $G$ and a single state $H$ (or some scaled versions of those)... unless you don't care about it being linear wrt $\alpha1, \alpha2,$ and $\theta,$ in which case they can merely be intermediaries in the time dependence of $A, B, C$ and $D$.

      Delete
    5. What I mean by "it will never be linear wrt $\alpha1, \alpha2$ and $\theta$" is that you will never be able to have a transfer function with those variables as inputs. Specifically, you'll never have anything like this: $H(s)/\alpha_1 (s)$.

      Delete
    6. Well, I guess if this was going to be easy, the Federal Reserve and central banks around the world would have figured it out by now. LOL

      I conceive this as a system with only one input $G$. Once $G$ starts, it continues until destroyed by my TR or SIM's $\theta $ . How how do we model a bump-in-the-road that persists over several time periods?

      The SIM model does that. If $G$ goes to zero but $\theta $ does not, over a few more time periods, all of the stored $G$ which we call $H$ wealth will disappear.

      It is almost like $G$ is the input and everything else is the output! It is all up-side-down, inside-out.

      BTW, the matrix code works just fine. It helps my thinking also. We have four matrices, each with as many points as terms in the longest equation. I knew that but still didn't have it internalized. Thanks!

      Delete
    7. "I conceive this as a system with only one input G."

      Good! That makes life easier.

      "How how do we model a bump-in-the-road that persists over several time periods?"

      I don't know what you mean by that. G can be anything as it stands.

      "It is almost like G is the input and everything else is the output! It is all up-side-down, inside-out."

      Well, except for the initial value of H, yes, that's correct. I don't see what's upside down or inside out about that though.

      "BTW, the matrix code works just fine"

      Do you mean A, B, C and D used as system matrices in Octave to model SIM?

      Delete
  4. Tom, I am bogged down at the moment. I see TWO delays in this state-space.

    1. G is in the state-space until end of period

    2. A part of G becomes H. H is also delayed in this state-space until end of period, creating the need for a second integrator. H will also need a feedback loop.

    I notice that $\theta $ and $\alpha 1$ are denominators. I think that makes them negative feedback so both would be degenerative.

    Following the same thinking, G and H are both numerators so would be positive feedback and regenerative. $\alpha 2$ however, is a fractional modifier of $H$ but, being fractional, it is also degenerative.

    I have not figured out how to display all that on a state-space representation.

    ReplyDelete
    Replies
    1. Roger, see my comment above. I changed my mind: I don't see a way to make this as a linear system with anything but the following: H (or a non-zero multiple of it) is your one and only state. Your state-space is fixed at one dimensional. And secondly, only G can be an input (or a non-zero multiple of it). The parameters, in contrast, can be time varying, but the result wrt any of them (or any combinations of them) is not a linear system: it will only be linear wrt G as an input and H as the state.

      The good news is that the problem is still easy to solve (as a time varying linear system). The bad news is that those parameters (or any functions of those parameters) cannot be treated as inputs to a linear system... there will never be any frequency domain analysis or transfer functions you can express them with. But you may not care about that.

      Delete
    2. ... but you can always write the system as a first order non-linear ordinary differential equation (ODE):
      $$\dot{x} = f(x,t)$$
      with $x$ a vector and some boundary condition $x(0) = x_0$.

      Delete
    3. This comment has been removed by the author.

      Delete
    4. In discrete time if you want a block diagram, just show $\alpha_1, \alpha_2$ and $\theta$ as inputs to $A_d$ and $C_d$ (I've added a subscript $d$ to signify discrete time and to distinguish the matrix $C_d$ from the output $C$), and $\alpha_1$ and $\theta$ as inputs to $B_d$ and $D_d$. Same goes for continuous time, except that $B_c$ (switching to a $c$ subscript for continuous time) is also dependent on $\alpha_2$. Even simpler is just to show them with a time dependence (since presumably $\alpha_1, \alpha_2$ and $\theta$ have a time dependence).

      Delete
  5. This comment has been removed by the author.

    ReplyDelete
  6. Roger, check it out: fun with MathJax: the discrete time and continuous time system matrices written as composite arrays:
    $$

    \left[

    \begin{array}{c|c}

    A_d & B_d \\

    \hline

    C_d & D_d

    \end{array}

    \right]_{SIM_d} =

    \left[

    \begin{array}{c|c}

    1-\frac{\alpha_2 \theta}{1 - \alpha_1 (1 - \theta)} & 1-\frac{\theta}{1 - \alpha_1 (1 - \theta)} \\[1.5ex]

    \hline

    \frac{\alpha_2}{1 - \alpha_1 (1 - \theta)} & \frac{1}{1 - \alpha_1 (1 - \theta)} \\[1.5ex]

    \frac{\alpha_2 \theta}{1 - \alpha_1 (1 - \theta)} & \frac{\theta}{1 - \alpha_1 (1 - \theta)} \\[1.5ex]

    \frac{\alpha_2 (1-\theta)}{1 - \alpha_1 (1 - \theta)} & \frac{1-\theta}{1 - \alpha_1 (1 - \theta)} \\[1.5ex]

    \frac{\alpha_2}{1 - \alpha_1 (1 - \theta)} & \frac{\alpha_1 (1-\theta)}{1 - \alpha_1 (1 - \theta)}

    \end{array}

    \right]

    $$
    $$

    \left[

    \begin{array}{c|c}

    A_c & B_c \\

    \hline

    C_c & D_c

    \end{array}

    \right]_{SIM_c} =

    \left[

    \begin{array}{c|c}

    \frac{1}{T_s} \log\left(1-\frac{\alpha_2 \theta}{1 - \alpha_1 (1 - \theta)}\right) & \frac{\theta + \alpha_1 (1-\theta) - 1}{T_s \alpha_2 \theta} \log\left(1-\frac{\alpha_2 \theta}{1 - \alpha_1 (1 - \theta)}\right) \\[1.5ex]

    \hline

    \frac{\alpha_2}{1 - \alpha_1 (1 - \theta)} & \frac{1}{1 - \alpha_1 (1 - \theta)} \\[1.5ex]

    \frac{\alpha_2 \theta}{1 - \alpha_1 (1 - \theta)} & \frac{\theta}{1 - \alpha_1 (1 - \theta)} \\[1.5ex]

    \frac{\alpha_2 (1-\theta)}{1 - \alpha_1 (1 - \theta)} & \frac{1-\theta}{1 - \alpha_1 (1 - \theta)} \\[1.5ex]

    \frac{\alpha_2}{1 - \alpha_1 (1 - \theta)} & \frac{\alpha_1 (1-\theta)}{1 - \alpha_1 (1 - \theta)}

    \end{array}

    \right]

    $$
    Where $T_s$ is the sample period of the discrete time model. Now can you write the 4x1 array of transfer functions in both cases (discrete and continuous time)?

    ReplyDelete
    Replies
    1. Perhaps a little less eye strain here.

      Delete
    2. ... no matter the number of states, inputs or output (i.e. the dimensions of x, u and y respectively) you can always write out the system matrices as one big rectangular matrix like that because their dimensions always conform to each other such that it's possible.

      Delete
    3. And where of course you can express the operation of the system with a matrix-vector multiply expression:
      $$

      \left[

      \begin{array}{c}

      x_{n+1} \\

      \hline

      y_n

      \end{array}

      \right] =

      \left[

      \begin{array}{c|c}

      A_d & B_d \\

      \hline

      C_d & D_d

      \end{array}

      \right]

      \left[

      \begin{array}{c}

      x_{n} \\

      \hline

      u_n

      \end{array}

      \right]

      $$

      $$

      \left[

      \begin{array}{c}

      \dot{x}(t) \\

      \hline

      y(t)

      \end{array}

      \right] =

      \left[

      \begin{array}{c|c}

      A_c & B_c \\

      \hline

      C_c & D_c

      \end{array}

      \right]

      \left[

      \begin{array}{c}

      x(t) \\

      \hline

      u(t)

      \end{array}

      \right]

      $$

      Delete
    4. I guess I should mention that in discrete time:
      $$x_n = H_n$$

      $$u_n = G_n$$

      $$y_n =

      \left[

      \begin{array}{c}

      Y_n \\

      T_n \\

      Y_{Dn} \\

      C_n

      \end{array}

      \right]

      $$
      and in continuous time:
      $$x(t) = h(t) $$

      $$u(t) = \dot{g}(t) $$

      $$y_m(t) =

      \left[

      \begin{array}{c}

      \dot{y}(t) \\

      \dot{t}(t) \\

      \dot{y_D}(t) \\

      \dot{c}(t)

      \end{array}

      \right]

      $$

      Delete
    5. where I added an $m$ subscript to the measurement vector $y_m(t)$ to distinguish it from the continuous time version of $Y$ (which I write as the rate $\dot{y}(t)$)

      Delete
  7. Nice Roger... you diagram is making progress.

    ReplyDelete
  8. This comment has been removed by the author.

    ReplyDelete
  9. Roger, I like your "One period time advance" block. If that's what I think it is i.e. a block that could be labeled like this:
    $\require{AMScd}$
    \begin{CD}
    x_{n+1} @<<< \bbox[yellow,border:2px solid black] {\underset{\text{by one}}{\overset{\text{Advance}}{Z^{+1}}}} @<<< x_{n}
    \end{CD}
    Instead of the usual
    $\require{AMScd}$
    \begin{CD}
    x_{n} @<<< \bbox[yellow,border:2px solid black] {\underset{\text{by one}}{\overset{\text{Delay}}{Z^{-1}}}} @<<< x_{n+1}
    \end{CD}
    Then I get what you're doing, and that's certainly a fair thing to do with a linear system conceptually. However it's not realizable... i.e. if we could do that in computers, they'd actually save time the more computations they did. ;)

    In the continuous time world, that's like an RLC circuit which starts to react before you apply an input voltage. Nothing wrong with it conceptually, but it is non-causual, and thus can't actually be constructed. Is that what you meant, or am I misreading you?

    ReplyDelete
    Replies
    1. This comment has been removed by the author.

      Delete
    2. This comment has been removed by the author.

      Delete
    3. This comment has been removed by the author.

      Delete
    4. This comment has been removed by the author.

      Delete
    5. Sorry for all the deleted comments Roger... looking at your diagram more, I think I got it wrong: you meant that block as a delay (as in my first diagram above), true? That way the output of the summing junction (on the right) being fed into the "advance" block would be $x_{n+1},$ but the output of the "advance" block (on the left) would be $x_n$, so it's really a delay block. That would make sense. I hope that's what you meant. Otherwise I'm confused.

      Delete
    6. "(as in my first diagram above)"

      should be

      "(as in my second diagram above)"

      just to be clear, this one:

      $\require{AMScd}$
      \begin{CD}
      x_{n} @<<< \bbox[yellow,border:2px solid black] {\underset{\text{by one}}{\overset{\text{Delay}}{Z^{-1}}}} @<<< x_{n+1}
      \end{CD}

      That's what your "advance block" is doing, right?

      (Shoot, I'd better go to bed... I'm making too many errors!)

      Delete
    7. Tom,

      Chuckles! I just got up for the second time and read quickly your comments! OK, a lot to say but first, I updated again.

      http://mechanicalmoney.blogspot.com/p/mathjax.html

      "On Delay Z-1 by one" , I do intend to advance the time. I wonder if Y is actually in Z-space? The state-space representation update speaks to that question, and points in that direction. (I am still trying to make every single piece fit into the puzzle.)

      We are not really predicting the future here. Instead, we are instantly updating the last best data to the logical present data.

      Can it be done? Yes, because we know the three steering factors (taxes, wealth tax, and rate of wealth draw down) which are all rates of change. We start with knowing the potential taxes which is "wealth". We convert wealth-and-G to Y (which must be a Z-transform but this is not solid in my mind)

      I wonder if the "wealth" conversion to Y causes an inversion in Z space?

      It will be a busy day at the ham fest for me. I travel to Yakima and will be back in the late afternoon. Will be thinking about something all the time!

      Delete
    8. This comment has been removed by the author.

      Delete
    9. Roger, here's something to think about: take a look at this diagram of a simple discrete time finite impulse response filter (FIR). We could replace all the $Z^{-1}$ blocks with $Z$ blocks if we reordered the taps, left to right, as $b_N, ..., b_0$ and we'd get the same output $y$ for the same input $x$, except with a time advance instead of a time delay. It works the same for linear system theory aside from the difference in timing... unfortunately it can't be built in reality. Imagine that the way it is now (with delays) that all the taps are 0 except $b_0 = 1$. So now you put in a signal $x$ and it immediately comes out as $y$ with no delay. But the filter with $Z$ blocks instead of $Z^{-1}$ blocks (and the taps in reverse order) would produce the output $y$ prior to (by $N$ samples) even applying the input $x$ because those $Z$ blocks are effectively each looking one time step into the future, so linked together in series like that, they can look $N$ samples into the future! ;)

      I can tell you, I've never in my life seen an actual $Z$ block, but if I could invent one I'd be rich!

      Normally with an equation like this:

      $$x_{n+1} = a x_n + b u_n$$

      The $x_{n+1}$ appears on the output of the summer. It's saying "this will be the value next time step." Then there's a delay register after that, so when the leading edge of the clock pulse arrives it "advances in time." I suppose it would be less confusing if they wrote the equation like this (which is exactly equivalent of course):

      $$x_{n} = a x_{n-1} + b u_{n-1}$$

      Then it's clear that the current value of $x$, called $x_n$, is dependent on past values of $x$ and $u$. It's true that the $Z$-transform of that equation is:

      $$z X(z) = a X(z) + b U(z)$$

      But that's of course equivalent to:

      $$X(z) = z^{-1} (a X(z) + b U(z))$$

      That equation says what's going on right there: a sum is formed, then it's delayed by one. Make sense?

      Of course with the feedback, what you're doing is more like an infinite impulse response (IIR) filter than an FIR filter, but still they are constructed with delays, not advances. Here's a 2nd order one for example.

      Delete
    10. You probably have seen these, but just in case, here's how a 1-bit wide $Z^{-1}$ (delay) block would actually be constructed. AFAIK nobody has yet invented a device to do the opposite!

      And just for completeness, here's one in the analog world (Ku band RF in this case): essentially just a long coil of shielded wire (or a wave guide in this case).

      Delete
    11. "It will be a busy day at the ham fest for me."

      Mmmm, ham. Sounds delicious! ;)

      Delete
    12. Tom, The ham fest was fun. I only spent $2.50 plus food and entry and gas and time. Some good visits with friends.

      Well, remembering that we can do things with math that we cannot do in the real world, I think my Eq. 10 is an

      "We could replace all the Z−1Z−1 blocks with ZZ blocks if we reordered the taps, left to right, as bN,...,b0bN,...,b0 and we'd get the same output yy for the same input xx, except with a time advance instead of a time delay. It works the same for linear system theory aside from the difference in timing... unfortunately it can't be built in reality.

      And yes, I think this is an amplifier with a mathematical advance. That makes the Y output an integral in Z-state-space which requires some explaining.

      First, notice that I updated Figure 1 yet again. This time I focused on the Y state-space. The text and vertical bars are intended to delineate the boundary of Y state-space existence. Y is a measurement that is unique to the period. Y is not something you can put your hands on, nor can you put Y into the bank. This makes me think that integral Y is the output from Z-space.

      Still more explaining but hold onto your hat. Wealth is just Y measured before wealth has been recaptured by tax. Looking at Y from a different way, if we doubled the time period, each increment of Y would be bigger but the increment to wealth would be smaller. Extending that analogy to the limit, Y would reach a level where each doubling of the time period would result in (near) zero Y increase and the wealth would be (near) zero.

      So, just to be clear, I think GDP is a measurement of the real world Z-state-space. Does that make any sense at all?

      I hope you held onto your hat!! LOL

      Next, I am going to try to get the same answer using matrix math. That may be a challenge for me due to my near zero matrix skills. I will try to fold in $\alpha 2$ which would produce a link not provided in my derivation of Equation 10.


      Delete
    13. Roger, I don't know what you mean by "Z-space."

      I put up another post myself on this... you may find it interesting. There's more to come, but I put some plots up for now. Basically I worked out in excruciating detail two things that have been bothering me:

      1. The canonical vs the non-canonical form. SIM is naturally in non-canonical form, but life is easier if you stay in canonical, and I find it often leads to confusion to drift away from canonical. Thus I work out the solution in both forms and compare what the variables look like side-by-side.

      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 bothered me for a long bit and is a bit more complicated to work out vs the canonical & non-canonical business. The canonical vs non-canonical business is no big deal, but it adds to the confusion with average vs instantaneous, so it helped to sort them both out at the same time.

      Here's the bit you may find interesting: I ran SIM at a MUCH higher sample rate (at about 1/100 the the sample period we usually use) with the continuous time equations. And I zoom in on those first few big time steps (corresponding to the discrete time steps) so you can see what's happening in gory detail.

      Here it is: enjoy!.

      Delete
    14. I only look at the the T output from SIM now... but I did work out the math for the rest. The plots are complex enough with just T though!

      Delete
    15. Roger, can you explain what you mean by "Z-space?" I don't understand you there.

      Delete
  10. Roger, I see you made edits to the figure, but I can't tell what they are. Did you replace it and then change it back?

    ReplyDelete
    Replies
    1. The edit to the picture was to change the text on the left side. Now Y -> T -> Yd, C, and H . I hope my 7:47 comment makes sense?

      Delete
  11. Tom, I added a Figure 0 and in Figure 1., on the far right side, changed
    Go to (Go-To). Just a few new lines of text at the top, above Figure 0.

    If you didn't understand my references to the Z plane, I am not using the analogy correctly. I am working on that. LOL

    ReplyDelete
  12. Tom, I added a Figure 0 and in Figure 1., on the far right side, changed
    Go to (Go-To). Just a few new lines of text at the top, above Figure 0.

    If you didn't understand my references to the Z plane, I am not using the analogy correctly. I am working on that. LOL

    ReplyDelete
    Replies
    1. Hi Roger, I see the figure, and yes, it makes sense. I think the jump you're referring to is due to the D matrix (AKA the "feed-through" matrix). If D = 0, then there'd be no jump in the measurements AKA outputs of the system.

      Hey, I thought of a real simple 1st order mechanical model that might help to understand. Imagine a system like the one pictured here with just a spring and a dashpot (dashpot = dampener, like a shock absorber). In this case both the spring and the dashpot have negligible mass. Now imagine an actuator connected to the right side of the spring which can compress or stretch the spring, and also that the left side of the dashpot is tied to an immovable wall, and that we're interested in the displacement of the canister side (right side) of the dashpot (call that displacement $x$) in response to disturbances from the actuator. Now say at time $t=0$ the actuator introduces a displacement of $u$ units of distance to the system. The force the spring produces is proportional to its compression from its resting quiescent state. We could write that as $F_s = k (u(t) - x(t))$ with $k$ a positive constant associated with the spring. The magnitude of the force the dashpot produces in resisting this movement is proportional to the rate of change (velocity) of point x wrt to the wall: $F_d = r v(t) = r \frac{dx(t)}{dt} = r \dot{x}(t)$ with $r$ a positive constant. And we know that these forces must always be equal and opposite, so we can write:
      $$r \dot{x}(t) = k (u(t) - x(t))$$
      or
      $$\dot{x}(t) = -\frac{k}{r} x(t) + \frac{k}{r} u(t)$$
      Substituting into our usual state space notation (where we have a single state here: the displacement $x$) we have:
      $$\dot{x}(t) = a x(t) + b u(t)$$
      Where $a = -k/r$ and $b = k/r$. If we can read the displacement $x$ directly from a horizontal measuring rod (say lying on the floor under the dashpot), then we can add a measurement equation (suppose the units of the measurement are equal to the units of $x$ and $u$):
      $$y(t) = x(t)$$
      and thus our measurement equation state space "matrices" (scalars in this case) are $c = 1$ and $d = 0$. How can we imagine a case where $d \ne 0$? Suppose for example that the actuator driving the spring on the right has coupled to it a second actuator that moves the measuring rod by a positive constant $q$ times the amount the actuator moves the spring. Then the measuring rod moves a distance $p = q u$ across the floor, and the dashpot displacement we read with this measuring rod is thus relative to the rod's movement: $y(t) = x(t) - p(t) = x(t) - q u(t)$. Now we have a non-zero $d = -q$ and our system is complete!
      $$\dot{x}(t) = a x(t) + b u(t)\\
      y(t) = c x(t) + d u(t)$$
      Our measurement $y(t)$ will behave exactly analogously to one of the measurements in SIM (i.e. $Y, T, Y_D$ or $C$). It will even have "jumps" in $y$ (cased by a non-zero $d$) in response to step changes in $u$. Of course it goes without saying that $u$ is analogous to $G$ and $x$ is analogous to $H$ (or rather their continuous time equivalents).

      What if you want a case where $a \ne -b$ or $c \ne 1$ or $d > 0$, well for the 1st of those I'd suggest looking at a slightly more complicated 1st order system with springs and dashpots such as this one (or alternatively reading the Wikipedia article associated with that picture and expressing the problem in terms stresses, etc, rather than displacements from an actuator like I did). If you want $c \ne 1$ then imagine the measuring rod is in different units than what $x, u$ and $p$ are measured in. And finally, if you want $d > 0$ then turn the measuring rod around $180^{\circ}$.

      Delete
    2. You might wonder "Why would anybody have a system with a coupled secondary actuator that moves the measuring rod?" ... Exactly!... that's one of the reasons $D$ is so often $0$ in real systems of interest. But I suppose there might be a known flaw in the measuring device that responds directly to the system's input in such a way (e.g. biasing the measurements) that introducing a non-zero $D$ might make sense. And then there's SIM of course! ;)

      Once you've mastered 1st order spring and dashpot systems (analogous to 1st order resistive-inductance (RL) or resistive-capacitive (RC) networks in basic circuit theory), then you can move on to 2nd order systems with a spring, a dashpot and a non-zero mass, like these. Notice how similar the discussion is there to the RLC circuit article that Jason linked to in one of his SFC related posts. Both systems are classic 2nd order systems, and in fact one of the lab exercises we used to have to do back in the day in electrical engineering classes was to use an analog computer to simulate a mechanical system with electrical components like RLC networks and Op Amps. (I'm going to guess that my class was probably the last one to use that in our lab... they were pretty antiquated at the time).

      And maybe you can think up more compelling practical situations in which $D \ne 0$. That one always makes me scratch my head.

      Delete
    3. This comment has been removed by the author.

      Delete
    4. ... (two comments up) for $d > 0$ it's probably better to imagine $q < 0$, i.e the 2nd coupled actuator... (perhaps nothing more than a hinged lever arm connected to the 1st actuator) moving the rod in the opposite direction to the spring. Otherwise turning the measuring rod around would also make $c = -1$. It'd be nice to be able to alter the value of $d$ without affecting anything else.

      Delete
    5. ... also notice that with both $k$ and $r$ positive constants, we get $a < 0$... meaning the system is stable, since the systems only pole will be at $s = a$ (and for continuous time systems, poles in the right half of the $S$-plane produce stable outputs. That's one way you can check yourself to make sure you wrote the equation correctly, because there's no way two passive elements (a spring and dashpot) are going to be unstable in this case. If you had some kind of magic spring that produced a force in the opposite direction from a normal spring when it was compressed or stretched, or a magic dashpot that produced a force aligned with velocity, then you could have an unstable system... and in fact neither would have to be "magic" ... you'd just need some sort of poorly designed (or rather poorly used) "smart" spring or dashpot that was active rather than passive.

      Delete
    6. How do I know the system's pole is at $s = a$? That goes back to the state update equation in our state space representation:
      $$\dot{x}(t) = a x(t) + b u(t)$$
      Taking the Laplace transform of both sides we have:
      $$s X(s) = a X(s) + b U(s)$$
      which means:
      $$
      \begin{align}
      X(s)(s - a) &= b U(s)\\[3ex]
      X(s) &= \frac{b}{s-a} U(s)\\[3ex]
      P(s) \triangleq \frac{X(s)}{U(s)} &= \frac{b}{s-a}
      \end{align}
      $$
      where $P(s)$ is the system's transfer function (from the input to the state anyway, but if $c=1$ and $d=0$, then from input to output as well). $P(s)$ obviously blows up when $s = a$, thus $a$ is the system's pole. $P$ stands for "plant" which is what control engineers call the given system under consideration, but $H(s)$ or $G(s)$ are often used as well (both of which would cause nothing but confusion in terms of SIM though!).

      Delete
    7. Let me know if my mechanical analogy makes sense. I figured you might like those best. In particular, did my explanation for how a jump could occur in the output due to an input actuator also moving the measuring rod make sense? It's not a very practical analogy, but I think it works. One good thing: all 1st order systems are alike! ;)

      Delete
    8. ... I'm confident you could come up with a better mechanical analogy than me: perhaps utilizing a mass and a dashpot for example (I'm pretty sure a mass and a spring won't work... that's getting into 2nd order territory).

      Delete
    9. ... a mass and a spring is analogous to an inductor and a capacitor (an LC circuit)... it's 2nd order and if there's really no dampening at all (no R or analogously no friction or dashpot), then it's a perfect oscillator as well, which will ring forever at it's characteristic frequency at constant amplitude without stopping (for it's impulse response and/or response to non-zero initial conditions (the same thing)). There are two poles in such a system (thus it's 2nd order), with equal magnitude weights occurring at conjugate locations on the imaginary axis at $\pm i \omega$ where $\omega$ is the characteristic frequency of oscillations in radians/sec.

      Delete
    10. Tom, This all makes a lot of sense to me. I followed it very well, but had to read it more than once.

      This about state-space representation, not seismometers!

      I can see why you find problems with my state-space Figure 1.

      Your main problem comes from the dx(t)/dy term where you expect a diminished response from the actuator due to the damping effects of the dashpot. I am claiming that the opposite is happening--we have an integral Sum dx(t)/dy.

      I also see a problem in my representation of D.

      I will puzzle about this longer. I think I am close to an explanation.

      Thank you very much.

      Delete
    11. Hi Roger,

      I'm glad that was helpful.

      I don't understand you here:

      "Your main problem comes from the dx(t)/dy term where you expect a diminished response from the actuator due to the damping effects of the dashpot. I am claiming that the opposite is happening--we have an integral Sum dx(t)/dy."

      But no matter: I'll let you stew on it a while and see what you come up with.

      Like I mentioned above if you have a "smart dashpot" or "active inverting dashpot" which produces a force vector in the opposite direction from a normal dumb passive dashpot, you'll destabilize this system: i.e. you'll have a > 0 and thus your pole will be in the right half plane. I still think the jump you're concerned about is entirely due to $D \neq 0$.

      Delete
    12. ... $D \neq 0$ just means a scaled version of the input (itself a scaled unit step in this case) shows up added into the output. That's why your outputs look like scaled version of this $(1 - e^{a t})$ where $a < 0$ added to scaled unit steps (or they're "sitting on a pedestal" so to speak). Think about what $f(t) = (1 - e^{a t})$ means with $a < 0$: it starts off at 0 and smoothly (not in a step!) climbs ever closer to 1 as t goes to $\infty$. That's what H is, and with $D = 0$ that's what all your outputs would look like (scaled version of it of course). But because $D \neq 0$ all those get added to constants (that's the "pedestal" that $f(t)$ sits on).

      Delete
    13. I see in the comments on my blog that you have a modified dashpot picture in your mind. Great, I will await your diagram! ;)

      Delete
    14. Also, I can't believe I wrote this above:

      "poles in the right half of the S-plane produce stable outputs."

      That should be "left half!!!"

      Delete
  13. ... also regarding "Z-space" that you mention above: remember the purpose of taking Fourier, Laplace or Z-transforms is to express a system in different coordinates, that we happen to call the "frequency domain." It's literally like using to set of geographic coordinates to express the location of items in your house: sometimes it's more convenient in one reference frame, and sometimes in another. So the system itself doesn't ever cross back and forth between reference frames... it doesn't care what reference frame we look at it in. That's up to us. But if we chose one frame (time domain or frequency domain), then it doesn't make sense to mix elements from both. We can cross back and forth, but you should never see something like this:
    $$\dot{x}(t) = A X(s) + B U(z)$$
    That's just horribly messed up!

    ReplyDelete
    Replies
    1. Let's get rid of that abomination:
      $$
      \require{enclose}
      \enclose{updiagonalstrike,downdiagonalstrike}{\dot{x}(t) = A X(s) + B U(z)}
      $$
      ;)

      Delete
  14. Roger, please add some more text to describe your dashpot figure 0.1 a little more. I don't think I understand it. Are there hinged sections in the middle?

    ReplyDelete
    Replies
    1. ... and are there two dashpots? One on the right and one on the left?

      Delete
  15. Roger, I've created a page with some questions for you here.

    ReplyDelete
  16. Also I started a new page with a diagram of a simple 1st order inductive-resistive (LR) circuit with just enough complexity to model SIM with 1 or 2 outputs. How are you at circuits? For me they're often easier to decipher than mechanical diagrams. Here it is if you're interested.

    The analogy I was thinking was the following:
    Vs is G
    i1 is H (the state)
    i2 and i3 are two distinct measurements (outputs).

    The diagram could be simpler and still give non-zero A, B, C and D state space constants, but in this case it gave just a enough flexibility to be interesting I thought.

    You want to give it a shot? Imagine Vs is a unit step of 20. I think I can write down the equations for this one, and pick values for R1, R2, R3 and L to match SIM and 1 or 2 outputs (and the state H as well). Like I say, it's a lot more clear for me, but I don't know how your circuit intuition is. You know that you cannot change the current through an inductor instantaneously, right? It all boils down to Kirchhoff's laws and this simple equation relating the current through and the voltage across an inductor:
    $$v = L \frac{di}{dt}$$

    ReplyDelete

Comments are welcomed but are moderated. It may take awhile before they appear to be viewed by all.