How does c2d(...,method) compute the matrices Ad and Bd for both 'zoh' and 'foh' as method?

Given a continuousTimeStateSpaceModel comprising the matrices A, B, C, and D (albeit D = 0) and given a sample time Ts, how does discreteTimeStateSpaceModel = c2d(continuousTimeStateSpaceModel,Ts,method) comprising the matrices Ad, Bd, C, and D work for both 'zoh' and 'foh' as method? As stated, I guess, C and D will remain unchanged in both cases, but how are Ad and Bd computed in each case?

i.e.

• 'zoh': Ad = ..., Bd = ...
• 'foh': Ad = ..., Bd = ...

Thanks in advance!

[Paul]

Hi Simon,

Can you provide a bit more insight as to the motivation for this question?

[Simon]

Hi Paul, sure. I use c2d(...,method) with both 'zoh' and 'foh' as method in the application of the theory of a thesis, e.g. for figures. However, in the text itself MATLAB never appears (No worries, I mention that I use it for the figures, ...) as the theory should be applicable with any programming language, thus I can't say that I used the c2d() command but want to specify the underlying formulas. But even if I indicated the c2d() command, I don't want to state that I used it without having a clue what it does internally. Hope this helps.

[Paul]

The Ad and Bd matrices in theory (assuming no delays in the model) for the ZOH case can be found in any number of places, e.g., bottom of page 9 at this link.

However, there is more than one option for numerically computing Ad and Bd for the ZOH case.

I suspect the situation is similar for the FOH case, but I'm not as familiar with that.

If you need to know the specific numerical method that c2d uses to compute Ad and Bd for state space models, the only way to know that for sure is to step through the code using the debugger for both cases of interest (assuming it's not hidden behind p-code).

[Simon]

I'm familiar with the analytical solution for Ad and Bd, and know of the existence of its various numerical equivalents, for the ZOH case. Similar to you, I have had the same assumption for the FOH case, so I guess, we have a similar knowledge level on that topic.

However, I couldn't find any documentation on how MATLAB performs the ZOH/FOH discretization using c2d().

Thanks to your hint about stepping through the code using the debugger, I (sort of) have narrowed down the discretization step to two functions called utDiscretizeZOH() and utDiscretizeFOH(), though it gets cryptic after that (at least for now and for me). So the next step would be to analyze what these functions actually do, while I would still prefer actual formulas over me analyzing the code and deriving formulas from it, as the code originates from formulas anyway. So any help is still appreciated.

[Paul]

IIRC, one approach for the ZOH case is as follows.

Let Mc = [A, B;Z1,Z2]*Ts

where Z1 and Z2 are zero matrices with number columns equal to that in A and B respectively and number of rows s.t. Mc is square.

Let Md = expm(Mc)

Then Md = [Ad, Bd;Z1,I2]

Example:

rng(100);
sysc = rss(3,1,2);
Ts = 0.1;
Mc = [sysc.a,sysc.b;zeros(2,3),zeros(2,2)]*Ts;
Md = expm(Mc);
Ad = Md(1:3,1:3);Bd = Md(1:3,4:5);
% compare
sysd = c2d(sysc,Ts,'zoh');
sysd.a - Ad

ans = 3×3

1.0e-15 * 0 -0.0555 0 0 -0.1110 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

sysd.b - Bd

ans = 3×2

1.0e-16 * -0.2776 0.0347 0 0.2776 0.0694 -0.0694

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

For the ZOH case, the C and D matrices would be unchanged. But that's not always true for the FOH case (it might depend on the actual structure of the plant). Here's an example where D changes from zero to non-zero

sysc = ss(zpk(-(1:2),-(3:5),1));
sysc.d

ans = 0

sysd = c2d(sysc,Ts,'foh');
sysd.d

ans = 0.0373

[Simon]

Instead of I2 you probably mean Z2.

And the respective difference's magnitude of 1.0e-15 and of 1.0e-16 should indicate that both approaches to the discretization practically produce the same result, right?

What you said about C and D in the ZOH/FOH case ist the one thing that I analyzed so far, as the two mentioned functions utDiscretizeZOH() and utDiscretizeFOH() include the following code (they call the continuous matrices a, b, c, d instead of A, B, C, D and the discrete matrices E, F, G, H instead of Ad, Bd, Cd, Dd):

• ZOH: E = ...; F = ...; G = c; H = d;
• FOH: E = ...; F = ...; G = c; H = ...;

[Paul]

No, I meant I2. The lower right, or the second (2), partition is an identity matrix (I).

rng(100);
sysc = rss(3,1,2);
Ts = 0.1;
Mc = [sysc.a,sysc.b;zeros(2,3),zeros(2,2)]*Ts;
Md = expm(Mc);
Md

Md = 5×5

1.2193 0.3476 0.2742 0.1780 -0.0296 -0.0208 0.7957 -0.0433 -0.0368 -0.1408 -0.4311 -0.5133 0.5139 -0.0270 0.0438 0 0 0 1.0000 0 0 0 0 0 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Simon]

Oh, okay, my mistake. Sorry!

Answer 1

[Pascal Gahinet]

Hi Simon

You are on the right track. ZOH and FOH is just analytic integration of the ODE:

dx/dt = A x(t) + B u(t)

under the assumption that u(t) is either piecewise constant (ZOH) or piecewise linear (FOH). In the ZOH case, u(t)=u[k] on the interval [k*Ts,(k+1)*Ts) so we have

d/dt [x;u] = [A B;0 0] * [x;u]

which integrates to

[x[k+1] ; u((k+1)*Ts)] = expm ([A B;0 0] * Ts) * [x[k] ; u[k]] .

We only use the first row since the values of u are given.

Things get more complicated with input and output delays and even more complicated with internal delays, but you got the idea.

Answer 2

[Simon]

Btw, I seem to have found out what c2d() does, yet, I'm not sure why it does what it does, at least in the 'foh' case.

First things first: Given a continuousTimeStateSpaceModel comprising the matrices A, B, C, and D (albeit D = 0) and given a sample time Ts, how does discreteTimeStateSpaceModel = c2d(continuousTimeStateSpaceModel,Ts,method) comprising the matrices Ad, Bd, Cd, and Dd work for both 'zoh' and 'foh' as method?

This can be answered as follows:

In general, the input vector u is part of the input vector space U and the output vector y is part of the output vector space Y:

'zoh':

'foh':

This was verified by comparing the results of c2d() and of the formulas above. This was done by computing their results and their difference in MATLAB, which indeed was exactly 0.

Now the tricky part where there are still some open questions:

How did I find out (again, I was looking for formulas and not for code, as the code concerned derives from formulas and not the other way round):

I consulted the documentation of c2d() and its discretization methods.

First, the 'zoh' case: The documentation doesn't include any information on what exact procedure is used in this case. Thanks to Paul's comments, I found this explanation (p.31-34). I guess you could argue to not include this in the documentation as this is trivial when working in and familiar with this field, however, I wouldn't suggest such claim/assumption.

Second, the 'foh' case: The documentation refers to a book that can be found online here. However, it refers to the wrong page, and when you eventually find the right one, the listed formulas describe a SISO (Single Input Single Output) instead of a more general MIMO (Multiple Input Multiple Output) case and/or are faulty, e.g. there are 0s where there should be 1s. This makes me wonder whether the referenced source really is the best one to explain the discretization approach taken. After all, assuming that the chosen discretization approach is valid, there likely exist multiple sources to refer to, that are correct throughout and explain it well.

To conclude, I have not yet understood entirely where the 'foh' discretization approach is coming from, including that I am missing an intuitve explanation why Dd wouldn't match D, especially when D was 0. I still hope for the clarification of the situation.

[Paul]

I believe the book you've linked to is Reference 2 on the linked doc page, but Reference 2 doesn't include a page number.

I'm looking at Reference 2, pages 205 and 206.

The listed formulas in eqns 6.42 - 6.47 are for SISO, but it seems like extending to MIMO should be straightforward, which looks like you've already done.

As for the "0's where there should be 1's" I guess your refering to the (2,3) term on the right hand side of eqn 6.44, which you're showing as an identity matrix?

After a quick look, I agree that the development from 6.42 to 6.47 could be clearer and include some more detail (and I'm questioning the location of the sampler in Figure 6.8).

[Simon]

Yes, the linked book is Reference 2. It includes a page number, yet not under References but in the text: "(see [2], p. 228)"

Why pages 205 and 206? I am especially looking at pages 160 and 161 (or 176 and 177, if you are referring to the pages of a pdf viewer).

For the 0's and 1's, yes, I am referring to the term you've described.

I'm glad you are saying that.

[Paul]

I was referring to a hardcover, hardcopy of the book. My third edition is dated 1998, not 1997 as in Reference 2 on the doc page. Not sure what you're seeing at that link. I'm not familiar with that site. Can that site legitimately provide full textbook downloads?

In any case, what you're looking at isn't what I'm looking at, and probably isn't the same as what the authors of the doc page were looking at, which likely explains the page discrepancies. Not sure how the third edition of a text book can be so different, but I'm not at all familiar with publishing practices.

[Simon]

Oh, I see. The link provided seems to be a backup for this site, which links to another site. They say that the full textbook can be downloaded, as the book went out of print, yet, they wanted to keep it available.

Okay, that at least explains the page discrepancies. Like you, I am wondering how one and the same edition of a textbook can be so different, however, no details are given regarding the referred to version of the book, e.g. hardcover, paperback, ebook.

[Paul]

What about the foh discretization isn't clear about what it's trying to do? Though I do think section 6.3.2 could be more expansive on the concept.

As far as the D-matrix goes, consider the following example to get an idea of what's happening

H = zpk([],[-2,-3,-4],1)

H =
 
          1
  -----------------
  (s+2) (s+3) (s+4)
 
Continuous-time zero/pole/gain model.

Following equation 6.39, we need to mulitply H(s) by 1/s^2 and then find the z-transform of the samples of the resulting impulse response

s = zpk('s');
T = 0.1; % for example
H1 = c2d(H/s/s,T,'imp')/T

H1 =
 
  3.4877e-06 z (z+8.301) (z+0.8353) (z+0.08405)
  ---------------------------------------------
    (z-1)^2 (z-0.8187) (z-0.7408) (z-0.6703)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

We see that even though H(s) had relative degree 3, H1(z) has a relative degree of 1. I think this situation is very common. I'm not sure if there are cases where H1(z) has a relative degeree more than 1. In any case, in this case the relative degree is 1.

The next step is to mutiply H1 by (z-1)^2/(Tz). Because that factor has two zeros one pole, it will reduce the relative degree of the prodcut to zero (numerator and denominator have same order)

z = zpk('z',T);
H2 = (z-1)^2/z/T*H1

H2 =
 
  3.4877e-05 z (z-1)^2 (z+0.8353) (z+8.301) (z+0.08405)
  -----------------------------------------------------
       z (z-1)^2 (z-0.8187) (z-0.7408) (z-0.6703)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

Cancel the common poles and zeros

H2 = minreal(H2)

H2 =
 
  3.4877e-05 (z+8.301) (z+0.8353) (z+0.08405)
  -------------------------------------------
       (z-0.8187) (z-0.7408) (z-0.6703)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

Verify it's the same result as c2d

c2d(H,T,'foh')

ans =
 
  3.4877e-05 (z+8.301) (z+0.8353) (z+0.08405)
  -------------------------------------------
       (z-0.8187) (z-0.7408) (z-0.6703)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

Because H2 has zero relative degree, its state space realization will have a non-zero D, even though the original continuous-time system did not have a direct feed through term

sys = ss(H2);
sys.D

ans = 3.4877e-05

For the zoh case, the multiplying factor is (z-1)/z, which will not change the relative degree of the system. So even though Z(H(s)/s) has relative degree 1

H1 = c2d(H/s,T,'imp')/T

H1 =
 
    0.00013338 z (z+2.988) (z+0.2134)
  --------------------------------------
  (z-1) (z-0.8187) (z-0.7408) (z-0.6703)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

The final system will still have relative degree 1 and therefore has a state space realization with D = 0.

Back to the foh case .... it will only have a realization with D = 0 if the degree of the numerator of H1 is at least two less then the degree of the denomintor of H1. Would be interesting to try to find a case, or the general conditions under which, where that characteristic holds true.

[Paul]

In looking into this question I also came across what I think is a peculiar behavior of c2d.

Define a transfer function as follows:

s = zpk('s');
Hc = 1/(s+1) + 1

Hc =
 
  (s+2)
  -----
  (s+1)
 
Continuous-time zero/pole/gain model.

Because Hc(s) is proper and has a direct feed path from input to output, I'm not sure the concept of impulse invariance makes sense because the continuous impulse response includes a Dirac delta function that can't be sampled. But c2d is happy to proceed

T = 0.1;
Hd = c2d(Hc,T,'imp')

Hd =
 
    0.1 z
  ----------
  (z-0.9048)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

and returns a result that is the discretization of the strictly proper portion of Hc(s).

c2d(1/(s+1),T,'imp')

ans =
 
    0.1 z
  ----------
  (z-0.9048)
 
Sample time: 0.1 seconds
Discrete-time zero/pole/gain model.

I'm surprised that c2d works this way instead of throwing an error, or at least a warning. As far as I can tell, this behavior is not documented.

[Simon]

Hi Paul,

I think I understand section 6.3.2 by now.

Nevertheless, I would have some comprehension questions or a series of them, that eventually build upon each other, if answering them individually is not to much to be asked for. If not, the following would be my first questions – I call them comprehension questions as I think I know their answers but want them to be verified:

How or why is the block diagram in Fig. 6.8 (its structure) deduced from ? To me, it looks like the objective was (and sort of poke at F(s)) to find a representation of F(s), that partitions F(s) into single operations that take meaningful entries and give meaningful results, as meaningful intermediate results. Is this structure of meaningful intermediate results prescribed by the structure of (6.42)? Is the structure of (6.42) chosen to form a system of first order linear differential equations (like the standard state-space representation)?

Like you, I am questioning the sampling in Fig. 6.8: Would the accurate structure of the block diagram be ?

Thanks in advance!

[Paul]

It seems to me that Figure 6.8 should be as you've described with the sampler being the first element as shown in Figure 6.5b.

I think that Figure 6.8 is shown that way to motivate eqn 6.42, as you said, which follows directly from the definition of the extrapoation filter, F(s). I'm not sure what you mean by "structure of (6.42) chosen." There really isn't much choice in how 6.42 is defined. Maybe 1/T term could be instead put on the integration from ubar to w, but that looks to be about it (and might complicate the rest of the derivation, which I haven't really dug into).

Have you tried looking for another source or two that has better exposition?

[Simon]

I probably could have clarified what I mean by "structure of (6.42)":

I've referred to the variable on the RHS of one ODE and its derivative on the LHS of the next ODE. However, I feel I understand your perspective of the deduction from Fig. 6.8 to (6.42), as it represents a better approach than mine which was just an attempt to reason Fig. 6.8. Therefore, overall, your answer corroborates my initial explanation for Fig. 6.8, nonetheless.

I've tried looking for alternative sources once, yet, I couldn't find something helpful immediately, so I focussed on the book again instead. Terms used for research would be "triangle-hold equivalent", "first-order-hold equivalent", "triangle approximation", and "ramp-invariant approximation" – they're mentioned either by the book or by the documentation.

Another question I have would be as follows:

Right before (6.45), the book states that the third equation of (6.42) can be integrated to show that . I doubt that, as I think it results , especially when I sketch or wdot in (6.42), respectively, and integrate it visually. If the book made a mistake here, I conceive a suspicion where that minus sign was accounted for. What would you get as a result for w[k]?

[Paul]

I didn't try to figure out what the authors meant by integrating "with care."

[Simon]

The general solution is

where .

I agree that

together with the FOH input

results

.

This resembles

.

(The book gets the sign in front of Gamma1 wrong, as it is not coherent with (6.45).)

This yields the results for v and w as shown in the book. However, shouldn't we get the same results when integrating wdot, as is the case for v? I would understand the book's deduction if it only was for the opposite sign, yet, I don't understand the time shift that's happening, as again, they define and I would define , at least from the integration.

[Paul]

I did not try to rederive all of the equations in 6.3.2, particuarly those for v[k] and w[k].

However, we can show by simulation that the equations in 6.3.2 appear to be correct, though there are several typos in 6.46 and 6.47, in addition to the one you found in the unnumbered equation just below 6.44.

Here we go ....

Define a harmless, second order, continuous-time system

A = diag([-1,-2]); B = [3;4]; C = [5,6]; D = 0;

Define the sample period and a 101-point time vector

Ts   = 0.1;
k    = 0:100;
tvec = k*Ts;

Define a random sequence of inputs

rng(100);
uvec = rand(size(tvec));

Define a function to evaluate u(t) with linear interpolation between the input sample points. With k*Ts <= t <= (k+1)*Ts, u(t) interpolates linearly between u[k] and u[k+1], i.e., our first-order, non-causal hold

u = @(t) interp1(tvec,uvec,t);

Use ode45 to integrate over the 10 discrete intervals with this FOH assumption

x(1,:) = [0 0];
for ii = 2:numel(k)
    [~,xtemp] = ode45(@(t,x) A*x + B*u(t),[tvec(ii-1),tvec(ii)],x(ii-1,:),odeset('MaxStep',Ts/100));
    x(ii,:) = xtemp(end,:);
end
y1 = (C*x.').';

Use lsim to simulate the system. Explicitly specify the FOH approximation

[y2,~,x] = lsim(ss(A,B,C,D),uvec,tvec,[0 0],'foh');

Now implement the equations from section 6.3.2. I derived these starting from 6.44 using the stated solution

v[k] = u[k]

and

w[k] = u[k+1] - u[k].

These match equation 6.47 in the text after correcting many typos in 6.46 and 6.47

F  = [A,B,zeros(2,1);0,0,0,1/Ts;0,0,0,0];
FT = expm(F*Ts);
Phi = FT(1:2,1:2); Gamma1 = FT(1:2,3); Gamma2 = FT(1:2,4);

Ak = Phi;
Bk = Gamma1 + Phi*Gamma2 - Gamma2;
Ck = C;
Dk = D + C*Gamma2;

Simulate these equations with lsim. It's important to understand that the state vector in these equations is NOT the state vector from the original system. Here, the state vector is z[k] = x[k] - Gamma2*u[k] (where x[k] is the state vector of the original system), so we need to adjust the initial condition in x[k] to the initial condition in z[k].

y3 = lsim(ss(Ak,Bk,Ck,Dk,Ts),uvec,tvec,[0 ; 0] - Gamma2*uvec(1));

Finally, simulate the c2d discretization, again accounting for the initial condition

y4 = lsim(c2d(ss(A,B,C,D),Ts,'foh'),uvec,tvec,[0 ; 0] - Gamma2*uvec(1));

All four results lie on top of each other

figure
plot(tvec,[y1,y2,y3,y4],'-o'),grid

Plot the differences. y2, y3, and y4 match exactly (at least to the precision of the plot), with the ode45 solution different by a very small amount.

figure
plot(tvec,y2-y3,tvec,y2-y4,tvec,y2-y1),grid

●

So it would appear that all three discrete-time approximations were all derived using w[k] = u[k+1] - u[k] as advertised, and they all very-well-approximate the solution to the continuous-time differential equation.

I did not try to see how these results for y3 would change if we derived the difference equations assuming w[k] = u[k] - u[k-1], but I don't think it will work out too nicely.