How does c2d(...,method) compute the matrices Ad and Bd for both 'zoh' and 'foh' as method?
Show older comments
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!
8 Comments
Paul
on 20 Jan 2025
Hi Simon,
Can you provide a bit more insight as to the motivation for this question?
Simon
on 20 Jan 2025
Paul
on 20 Jan 2025
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).
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
sysd.b - Bd
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
sysd = c2d(sysc,Ts,'foh');
sysd.d
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
Simon
on 21 Jan 2025
Answers (2)
Pascal Gahinet
on 29 Jan 2025
1 vote
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.
Simon
on 29 Jan 2025
0 votes
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):
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.
12 Comments
Paul
on 29 Jan 2025
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
on 30 Jan 2025
Paul
on 30 Jan 2025
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
on 30 Jan 2025
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)
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
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
Cancel the common poles and zeros
H2 = minreal(H2)
Verify it's the same result as c2d
c2d(H,T,'foh')
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
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
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.
Define a transfer function as follows:
s = zpk('s');
Hc = 1/(s+1) + 1
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')
and returns a result that is the discretization of the strictly proper portion of Hc(s).
c2d(1/(s+1),T,'imp')
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
on 3 Feb 2025
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)?
? 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!
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?
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]?
. 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
on 4 Feb 2025
I didn't try to figure out what the authors meant by integrating "with care."
Simon
on 4 Feb 2025
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.
and I would define , at least from the integration.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);
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.').';
[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.
Categories
Find more on Aerospace Applications in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
